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PREFACE 


This report describes the status of work being performed on a Novel 
Approach in Formulation of Special Transition Elements research, (Grant NAG 3- 
790). The research is being monitored by Dr. C.C. Chamis of the NASA Lewis 
Research Center. The graduate research assistants for this term were Mr. 0. 
Odabas and Mr. M. Yahiaoui. Mr. Odabas and Mr. Yahiaoui, for whom I am the 
dissertation advisor, are both Ph.D. degree students in the Aeronautical and 
Astronautical Engineering Department. 


Nesrin Sarigul 



I. INTRODUCTION 


The objective of this research program is in the development of more 
accurate and efficient advanced methods for solution of singular problems encoun- 
tered in various branches of mechanics. The research program can be categorized 
under three levels. First two levels involve with the formulation of new class 
of elements called "Mesh Interface Elements" (MIE) to connect meshes of tradi- 
tional elements either in three dimensions or in three and two dimensions. The 
finite element formulations are based on the boolean sum and blending operators. 
This report describes the results obtained from the first two levels of the pro- 
gram. It may be noted that, at present, the second level of the program is 
being conducted under NAG 3-790. 

In today's advanced aircraft and space structure applications, steep tem- 
perature and/or stress gradients are commonly encountered. The analysis methods 
need to incorporate these steep gradients into the solution efficiently and 
accurately. Mesh Interface Elements are being formulated and tested in this 
research to account for the steep gradient effects. At present, the heat trans- 
fer and structural analysis problems are being formulated from uncoupled theory 
point of view. 

The status report, first, summarizes the general formulation for heat 
transfer and structural analysis by including the newly introduced varying 
material properties at material nodal points of the elements concept. Then the 
formulation of mesh interface elements are detailed. On the computational 
efficiency side, a hidden- symbolic computation concept developed by the author 
is given. Verification examples are included from heat transfer and structural 
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analysis problems. Appendix includes listings of the computer modules that are 
developed for this purpose. 
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II. FINITE ELEMENT FORMULATION FOR HEAT TRANSFER AND STRUCTURAL ANALYSIS 

Thermal effects induced by aerodynamic heating on advanced aircraft and 
spacecraft systems of current technology requires special analysis procedures 
in order to design these structural components to fulfill the specific mission 
requirements. The steep temperature and/or stress gradients and unusual advanced 
geometry and material concepts are being major items to deal with. In this 
research , the steep gradients problems are tackled by utilizing the Mesh Inter* 
face Elements and variable material properties at material nodal points elements 
concept . 

II. a. Heat Transfer Analysis 

The heat conduction equation in three dimension is 

pc p 21 “ _1 [k*, SI + SI + Kz 81 ] + L. [*yx -3_+ kyy £T + ky 2 il] 

at ax ax ay az ay ax ay az 

+ [k„ dT + k iy £T + k 22 21 ] + f 

dz 3x dy dz 

where p is the mass density of the material, Cp is the heat capacity, T is the 
absolute temperature, k^, . . .k 22 are the heat conduction coefficients and f is the 
heat source term. 

The equation will be solved subjected to the boundary conditions in terms 
of prescribed either the temperature and/or temperature gradients. 

For homogeneous anisotropic solids the equation becomes 

PC. 51 - k„ Sh. + kyy dll + k 22 + (k^y + ky,) ^T_ ( ky , + ^y ) 

dt ax 2 ay 2 az 2 axay ayaz 


+ (k^ + k„ £fl_) + f 

azax 
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in case of isotropic material we have 


pc T - k (2fl + + £t ) + f 

ax 2 ay 2 a z 2 

for the steady- state conditions, we obtain the Poisson's equation as 
k 2 T + f - 0 

where 

v 2 - ( d 2 + d 2 + a 2 ) 
ax 2 ay 2 az 2 

By utilizing the method of weighted residual on the heat conduction equation, 
we obtain 

- JpCpW dV + fwk( . T) dV + JwfdV - 0 
v 3T ^ ^ 

where W is the weighting function and T is the approximated temperature field 
and V is the volume of the structure. By utilizing the Green-Gauss theorem on 
the first term of the equation above, we obtain 

Jk. W. T dV - Jw.fdV + fn.k. TdS 

V 4 V 4 S 

where S is the surface of the domain. 

The discretized equations for NE number of subdomains may be given as 

HZ f HZ r M , 

2 Jk. W. TdV - 2 fw.f.dV. + 2 fn. TdS 

where V # and S t are the element volume and surface mean, respectively. 

By utilizing the shape functions used in approximating the temperature 
field to be the weighting functions, we achieve Galerlin type formulation. At 
an element level the equations are 
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[K C ]M)* - (Q)* + {Q s l 


where [K e ]* is the finite element conduction matrix, { Q ) * is the finite element 
load vector and {Q,}* is the heat source vector. 

It may be noted that material properties are assumed known at material 
nodal points within an element as it is detailed in the following section. 

The present numerical applications deal with the steady state equations 
and isotropic material behavior for the time being. 

II. b. Structural Analysis 

The principle of virtual work will be employed, 

|i| 

J(5e) T odV - J(5A) T bdV + J(SA) T sds + E^ (^A) 1 ?, 

where (fie) T is the virtual strain vector, (SA) T is the virtual displacement 
vector, b is the body force vector and s is the prescribed surface tractions 
vector and ? L indicate the point loads. 

For an anisotropic material by utilizing the generalized Hooke's law the 
equation above can be written as 

J(Se) 1 . A.e dV - J(5e) T .A.e. init dV - J(fiG) T o lnlt dV +J(SA) T .bdV + J\$A) T .s.dS 
+ E (fiA) 1 ?! 

i«l 

However, by discretizing the region to be analyzed into an NE number of 
finite elements, we obtain for one element (e) 

f(«€) T .A.edV - f(Se) T . A. e init dV - J(«e) T o lnlt dV + f(«A) T .b.dV + J(SA) T .s.dS 
J v« v« v« s« 

+ (E (fiA^P,)* 

t-i 
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where V. and S. denote the volume and the boundary surface of the element. 

The field variable A, displacement field, will be approximated in terms 
of the nodal displacements times the shape functions. In addition, the material 
properties are assumed to be known at "material nodal points," NM as introduced 
in this research, hence, we have 

A - Z NMATi . A i 

where NMAT i are the "material shape functions." It may be noted that it is 
possible to use material shape functions to be the same as the field variable 
shape functions. This approach, newly introduced here, will increase solution 
efficiency for problems involving large gradients. 

A summary of the equations takes the form of 

M.A + K. A - P 

It may be noted that unlike the traditional formulations, in calculation of the 
stiffness and mass matrix material properties are considered to be known at 
material nodal points. Therefore, the stiffness matrix, for example, will be 
calculated as 

m 

B t (E NMAT^A*) .B dV 

VC i-1 

Application problems include three-dimensional mesh interface elements 
under different loading conditions. The results obtained from test examples 
ar i reported in numerical examples section. 


K«- J 
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III. FORMULATION OF MESH INTERFACE ELEMENTS 


III. a. General Concepts on Element Formulation 

Formulation of the "Mesh Interface Elements" (MIE) are based on the boolean 
cum in dealing with analysis aspects in three dimensions, it is known that the 
higher order elements are too costly to use throughout the domain. It may be 
noted that, we need these higher order elements only at certain local areas of 
the domain. The rest of the domain can be modeled by utilizing lower order 
elements. In this research, a series of Mesh Interface Elements are being 
developed to connect meshes of different types of elements. 

The boolean summation in three dimensions is given as 

P[F] - P r ©(P,©P t ) [F] 

where P r , P s and P t are* the projector operators, and F is the field variable. 
The boolean summation yields with an approximation on the field variable F so 
that along the boundaries of the domain the continuity requirements are met in 
an exact manner. 

IH.b. Three-Dimensional Mesh Interface Elements 

In order to demonstrate the power of new Mesh Interface Elements, the most 
commonly used three-dimensional element meshes are considered. These are the 
meshes of 8-noded and 20-noded elements. A three-dimensional Mesh Interface 
Element then is formulated by utilizing the boolean sum and the projector 
operators . 
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Figure 1. Three-Dimensional 12-noded Mesh Interface Element 

(12-MIE) 

The formulation domain is [0,1] x [0,1] x [0,1]. This element will connect 
to an 8-noded solid element from (r,o,t) surface and to a 20-noded solid element 
from (r.l.t) surface. In addition, the field variable for this research, either 
tbe temperature field or the displacement field, is required to be continuous 
on the boundary surfaces of the element. 

In order to construct the element to meet the constraints above, the 
projector operators become, 

P r [F] - rF(l,s,t) + (l-r)F(o , s , t) 

P.[F] - sF(r , 1 , t) + (l-s)F(r ,o, t) 

P t [F] - tF(r , s , 1) + (l-t)F(r,s,o) 



The boolean sum then takes the form of 

P[F] - (P r + P. + P t - P r P, - P.P t * Vt + P r P,P t )[ p ] 

Now, by selecting the surface interpolation functions , F(l,s,t)... F(r,s,o) 
to satisfy the continuity requirements we obtain the approximation function as 

F “ P - P[F] 

or 

P - (r-rs) (l-t)F 1 + (2rst 2 - 2r z st + 2r 2 s - rst - rs)F 2 
+ (-2rst 2 - 2r z st + 2r z s + 2st 2 + 5rst - 3st - 3rs + s)F 3 
+ (1-s-r+rs) (l-t)F 7 + (rt-rst)F 5 + (2rst 2 + 2sr 2 t - 3rst)F 6 
+ (2r z st - 2rst 2 + 2st 2 - rst - st)F 7 + (t-rt) (l-s)F 8 
+ (4rs-4sr 2 ) (1- t)F 9 + (4rst-4rst z )F 10 + (s-rs) (4+-4t z )F u 
+ (4rst-4r 2 st)F 12 

where F lf ...F 12 are n °dal values of the field variable. 

The material properties are considered to be known at nodal points of the 
element and 

NM 

[A] - E NMAT^ 
i-1 

where NMATi are the material shape functions. 

The stiffness matrix is obtained and computer modules are included in the 
Appendix . 

This element is tested for both heat transfer and structural analysis 
applications and performs well as a Mesh Interface Element. 
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III.b.2 Three -Two Dimensional Mesh Interface Elements 

These elements are at present being tested and the results, together with 
the detailed formulation, will be reported at the end of the second level 
research. 
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IV . COMPUTER MODULES 


Mesh Interface Elements are coded to solve heat transfer and structural 

/ 

analysis problems. Listing of these computer modules are included in the 

Appendix, 

IV. a. Description of Computer Modules 

ELEM Constructs the element stiffness matrix for 8-, 12- and 20-noded 
three-dimensional elements. 

ELEM08 Calculates the derivatives of the shape functions of the 8-noded 
element . 

ELEM 12 Calculates the derivatives of the shape functions of the 12-noded 
mesh interface element. 

ELLM20 Calculates the derivatives of the shape functions of the 20-noded 
three-dimensional element. 

SHAP08 Calculates the shape functions of the 8-noded element. 

SHAP12 Calculates the shape functions of the 12-noded mesh interface 
element. 

SHAP20 Calculates the shape functions of the 20-noded element. 



Figure 2. Element Numbering Scheme For Three-Dimensional Elements 
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In addition to the modules given on the previous page, a general purpose 
assembly module is structured in order to include two- and three-dimensional 
elements together with the mesh interface elements. Further, a "hidden symbolic 
computation" scheme is developed in order to compute the integrals, resulting 
from the element formulations, exactly. A library of these modules are: 

POLDIF Differentiates a given polynomial 
POLMLT Multiplies two polynomials 
POLADD Adds two polynomials 
POLINT Integrates a given polynomial 

POLIEV Evaluates an integral at its upper and lower limits 

It may be noted that the hidden symbolic computation scheme works well 
except that it is a bit time consuming at this stage to convert everything in 
this form. 

At present, due to the time limitations, and also due to the maj^or 
objectives of the research being rather different from the hidden- symbolic 
computation concept, this approach is set aside for the time being. 

IV. b. List of Input Variables 

NNODES: Number of nodes 

NELEM: Number of elements 

NBOUND: Number of boundary surfaces 

COOR (I , J) : x,y,z coordinates of the Ith node 

NC0NE(I,1): Element type for the Ith element 

(NCONEC( I , J ) , J-2 , NC0NEC(I,1)+1): Nodal connectivities 


12 


FORCE (I): Element forces 


NBKIND(I,1), NBKIND(I,2), NBKIND(I,3): Surface number, number of nodes, 
element number for the Ith boundary surface 



Figure 3. Surface Number for Three-Dimensional Elements 

CBOUN (1,1), CBOUN (1,2), CBOUN (1,3): Prescribed heat flux convection 

coefficient outside temperature 

COEF: Conduction coefficient 

NTEMP: Number of nodes with prescribed temperatures 

ITEMP, TEMP(ITEMP): Node number, corresponding temperature 

IV. c. Computational Considerations 

The computations are being performed on a VAX 11/780 computer. Larger 
scale applications will be solved at NASA Lewis Cray-XMP. It may be noted that 
since Cray uses the VAX 11/780 type computers as a front end, Cray runs will be 
performed readily. 
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V. Numerical Verifications and Concluding Remarks 


Mesh Interface Elements are being tested and the results obtained from 
three-dimensional mesh interface elements are included in this section. 
Numerical examples include heat transfer and structural analysis problems by 
using different mesh sizes and different element types. It may be noted that 
very good results are obtained when a mesh interface element is used to connect 
different element meshes together. 

In concluding, the research will continue as it is given in the proposal. 
On the computational side, at this time, we need to utilize the Cray computer 
at NASA Lewis Research Center. 
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V.a. Heat Transfer Applications 
STEADY STATE HEAT CONDUCTION 3 ELEMENT SOLUTION 


A ■*- 


\ 



a 2 T + afr + a*T - o 

3x 2 3y 2 dz 2 


Boundary Conditions 


T|jc-o 

- y 2 - 3z 2 

fill -4 



3 x li-i 

H 

1 

o 

- 2x 2 - 3z 2 

fill - 6 


3y|y_ 3 

T|*-o 

- 2x 2 + y 2 

5T| - -6 


I z-l 
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Mesh III. 12-noded/12-noded/12-noded elements 
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FINITE ELEMENT MESHES (Cont.) 



Results from each mesh types given in terms of temperature distributions 
and also % error as compared to the analytical values . Mesh interface element 
performs well as it is seen from the nodal temperatures. 
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** TEMPERATURE DISTRIBUTION ** 


WV6-SW X C »/ 12./20*) 


NODE NUMBER 


TEMPERATURES 


1 

0 . 2000G00000D+0 1 


0.72101 19546D— 28 

3 

1000000000D+01 

4 

- . 3000000000D+0 1 

5 

0. 3000000000D+01 

6 

0. 1 OOOGOOOOOD+O 1 

7 

0.64B9107591D-01 

8 

20000 OOOOOD+O 1 

9 

0 . 5983739838D+0 1 

10 

0. 4500C OOOOOD+O 1 

11 

0. 4000000000D+01 

12 

0. 5378099O78D+O 1 

13 

0 . 32500 OOOOOD+O 1 

14 

0 . 2907929283D+0 1 

IS 

0. 13085231 16D+01 

16 

0. 1000000000D+01 

17 

0 . B250000000D+0 1 

IS 

0. 6250000000D+01 

19 

0. 52248367950+01 

20 

0. 32500000000+01 

21 

0. 1 100000000D+02 

*> 

*»<(» 

0. 9500000000D+0 1 

23 

0. 9000000000D+0 1 

24 

0. 10181 67 1 75D+02 

4 bW 

0. B250000000D+01 

26 

0. 80 10500 1060+01 

27 

0 . 655666 1 842D+0 1 

28 

0 . 6000000000D+0 1 
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f.} i ' H ti rJ f.* r ) rj r.» 


NODE NUMBER 


iUMERICAL VALUES 


ANALYTICAL VALUES 


ERROR 


4 

* 

V.200000000D+01 

0 . 200000000D+0 1 

0.0000 7. 

4m 

0. 72101 1955D-28 

0. OOOOOOOOOD+OO 

0.0000 *** 

3 

i ooooooood+o i 

— 0. 1 OOOOOOOOD+O 1 

0.0000 7. 

4 

-C» . 300000000D+0 1 

— 0. 300000000D+0 1 

0.0000 7. 

5 

0.300000000D+01 

0 . 300000000D+0 1 

0.0000 7. 

<3 

•:<. i ooooooood+o 1 

0 . 1 OOOOOOOOD+O 1 

0.0000 7. 

/ 

0. 64B910759D— 01 

0 . OOOOOOOOOD+OO 

0.0649 *** 

s 

-0. 200000000D+01 

-0 . 2000000000+0 1 

0.0000 7. 

o 

0.598373984D+01 

0 . 600000000 D+O 1 

-0.2710 7. 

10 

0.450000000D+01 

0. 450000000 D+O 1 

0.0000 7. 

* « 
m • 

0. 400000000D+01 

0 . 400000000 D+O 1 

0.0000 7. 

22 

0 . 537809908D+0 1 

0. 525000000D+01 

2.4400 7. 

13 

o. 325000000 D+O 1 

0. 323000000D+0 1 

* 0.0000 7. 

1 4 

C* . 290792928D+0 1 

0.300000000D+01 

-3.0690 7. 

IS 

0. 1 308523 12D+01 

0. 1 50000000D+0 1 

-12.7651 7. 

16 

0. 1 OOOOOOOOD+O 1 

0. 1 OOOOOOOOD+O 1 

0.0000 7. 

17 

0. 825000000D+01 

0 . 825000000D+0 1 

0.0000 7. 

IS 

0.625000000B+01 

0. 625000000D+01 

0.0000 7. 

1? 

0.522483679D+01 

0 . 525000000D+0 1 

-0.4793 7. 

20 

0 . 325000000D+01 

0. 325000000D+01 

0.0000 7. 

4 

*> 4 

C*. 1 1 OOOOOOOD+02 

0. 1 1 OOOOOOOD+02 

0.0000 7. 


0.930000000B+01 

0.950000000D+01 

0.0000 7. 

23 

C*. 900000000D+01 

0. 900000000D+01 

0.0000 7. 

24 

0. 101816718D+02 

0. 102500000D+02 

-0.6666 7. 

*»w 

0. 825000000D+01 

0 . 825000000 D+O 1 

0.0000 7. 

2* 

0. 80105001 1D+01 

0 . 800000000D+0 1 

0. 1313 7. 

#* *^i 

*» / 

0. 655666 1B4C+01 

0.650000000D+01 

0.8717 7. 

23 

0.600000000D+01 

0 . 600000000D+0 1 

0.0000 7. 


M6S» X 
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♦ * TEMPERATURE DISTRIBUTION ** AL 


NODE NUMBER 


TEMPERATURES 


1 

0. 20000 00 OOOD+O 1 

2 

0. 3000000000D+00 

3 

0. OOOOCOOOOOD+OO 

4 

0. 1 230000 000D+0 1 

3 

7500000000D+00 

6 

1000000000D+01 

7 

- . 2300000000D+0 1 

8 

- . 3000000000D+0 1 

9 

0. 22300 OOOOOD+Ol 

10 

0. 23000 OOOOOD+OO 

22 

73000000000+00 

12 

- . 2 730000 OOOD+O 1 

13 

0 . 30000000000+0 1 

14 

0. 1 300000000D+0 1 

13 

0. 1000000000D+01 

16 

0 . 22300 OOOOOD+O 1 

17 

0. 2300000000D+00 

18 

0 . 3866783466D- 1 4 

19 

1 SOOOOOOOOD+O 1 

20 

- . 2000000000D+0 1 

21 

0. 4250000000D+0 1 

• ^ 

0 . 22300000000+0 1 

23 

0. 1 230000000D+0 1 

24 

750 OC OOOOOD+OO 

23 

0 . 6000000000D+0 1 

26’ 

0 . 45000000000+0 1 

27 

0. 4000000000D+0 1 

23 

0 . 3250C OOOOOD+O 1 

29 

0. 32300000000+01 

30 

0 . 3000C00000D+0 1 

31 

0. 1 500000000D+0 1 

32 

0. 1000000000D+01 

33 

0 . 82500 OOOOOD+O 1 

34 

0. 6230000000D+01 

33 

0. 3250C00000D+0 1 

36 

0. 3230000000D+0 1 

37 

0. 1 l 00C00000D+02 

39 

0. 9500000 OOOD+O 1 

39 

0. 9000000000D+0 1 

40 

0. 1 023C00000D+02 

41 

0. 8230000000D+0 1 

42 

0. 8000000000D+0 1 

43 

0. 6300000000D+0 1 

44 

0. 6000000000D+0 1 


Zo /2- /Zo'i 
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NODE NUMBER NUMERICAL VALUES ANALYTICAL VALUES ERROR 


« 

■:>. 200000000D+01 

0. 200000000D+01 

0.0000 7. 

** 

0. SOOOOOOOC'D+OO 

0. 500000000D+00 

0.0000 7. 

3 

O.OOOOOOOOOD+OO 

O.OOOOOOOOOD+OO 

0.0000 *+* 

4 

0. 125000000B+01 

0. 125000000D+01 

0.0000 7. 

ar 

-0 . 730000000B+00 

-0. 750000000D+00 

0.0000 7. 

6 

-0. 1 OOOOOOOOD+O 1 

-0. 1 OOOOOOOOD+O 1 

0.0000 7. 

/ 

-0 . 250000000D+0 1 

-0.250000000D+01 

0.0000 7. 

8 

-0. 300000000D+01 

-0. 300000000D+0 1 

0.0000 7. 

? 

0 . 225000000D+0 1 

0.225000000D+01 

0.0000 7. 

10 

0.25000000C-D+00 

0.250000000D+00 

0.0000 7. 

• • 
A • 

-0 . 750000000D+00 

-0. 750000000D+00 

0.0000 7. 

« * 

-0 . 275000000D+01 

-0. 275000000D+01 

0.0000 7. 

13 

0 . 300000000D+0 1 

0 . 300000000D+0 1 

0.0000 7. 

34 

0. 1 50000000D+0 1 

0. 1 50000000D+0 1 

0.0000 7. 

15 

0. 1 OOOOOOOOD+O 1 

0. 1 OOOOOOOOD+O 1 

0.0000 7. 

1* 

C 1 . 225000000D+01 

0 . 225000000D+0 1 

0.0000 7. 

17 

0. 250000000D+00 

0 . 250000000D+00 

0.0000 7. 

IS 

0 . 58667 8547D— 1 4 

O.OOOOOOOOOD+OO 

0.0000 *** 

1? 

-C>. 1 50000000D+0 1 

-0. 1 50000000D+0 1 

0 . 0000 /. 

20 

200000000D+01 

-0. 200000000 D+0 1 

0.0000 7. 

« 

*• — 

0.425000000D+01 

0 . 425000000D+0 1 

0.0000 7. 

~ •» 

0 . 225000000D+0 1 

0 . 225000000 D+0 1 

0.0000 7. 

A.O 

0. 125000000D+01 

0. 1 25000000D+0 1 

0.0000 7. 

24 

-0 . 750000000D+00 

-0. 750000000D+00 

0.0000 7. 

~*ar 

0.600000000D+01 

0 . 600000000D+0 1 

0.0000 7. 

id 

0.450000000D+01 

0 . 450000000D+0 1 

0.0000 7. 

- / 

0 . 400000000B+0 1 

0 . 400000000D+0 1 

0.0000 7. 

23 

0.525000000B+01 

0 . 525000000D+0 1 

0.0000 7. 

29 

0.325000000B+01 

0 . 325000000D+0 1 

0.0000 7. 

30 

0.300000000D+01 

0. 300000000D+01 

0.0000 7. 

3 1 

C=. 150000000D+01 

0. 150000000D+01 

0.0000 7. 

32 

0. 1 OOOOOOOOD+O 1 

0.1 OOOOOOOOD+O 1 

0.0000 7. 

23 

0 . 825000000B+0 1 

0 . 825000000D+0 1 

0.0000 7. 

34 

0.6250000002+01 

0 . 625000000D+0 1 

0.0000 7. 

•«r — # 

0.525000000D+01 

0 . 525000000D+0 1 

0.0000 7. 


0 . 325000000D+0 1 

0 . 325000000D+0 1 

0.0000 7. 

. •% * 

0. 110000000D+02 

0. 110000000D+02 

0.0000 7. 

33 

0 . 950000000D+0 1 

0. 930000000D+01 

0.0000 7. 

29 

•0 . 900000000B+0 1 

0 . 900000000D+0 1 

0.0000 7. 

*0 

0. 1 02500000D+02 

0. 1 02500000D+02 

0.0000 7. 

41 

0 . 825000000D+0 1 

0 . 825000000D+0 1 

0.0000 7. 

*2 

0 . 800000000D+0 1 

0. 800000000D+01 

0.0000 7. 

h3 

0. 650000000D+01 

0 . 650000000D+0 1 

0.0000 7. 

1A 

0 . 600000000D+0 1 

0. 600000000D+01 

0.0000 7. 


MbdV* H Cz®/2 •/* 
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** TEMPERATURE DISTRIBUTION ** HC (v'l./ix Ax') 


NODE NUMBER TEMPERATURES 


1 

0. 2000000000D+0 1 

2 

0.46721 71 732D-28 

3 

1000000000D+01 

4 

- . 3000000 OOOD+O l 

3 

0. 3000000000D+01 

6 

0. 1 500000000D+0 1 

7 

0. 1000000000D+01 

8 

0. 249754643 1D+01 

9 

0. 23000000000+00 

10 

-. 17049343770+00 

11 

-. 1849798702D+01 

12 

2000000000D+0 1 

13 

0. 6000000000D+0 1 

14 

0. 4000000000D+0 1 

IS 

0.31 37206428D+0 1 

16 

0. 1000000000D+01 

17 

0. 1 100000000D+02 

18 

O '. 9500000000D+0 1 

19 

0. 9000000000D+01 

20 

0. 1049921 136D+02 

21 

0 . 8230000000D+0 1 


0. 7822844036D+0 1 

23 

0.61318664230+01 

24 

0. 60000000000+0 l 


22 


NUMBER 

NUMERICAL VALUES 

ANALYTICAL VALUES 

ERROR 

* 

* 

0 . 20000000C D+0 1 

0 . 200000000D+0 1 

0.0000 

*■> 

0.46721 71 73B-28 

0. OOOOOOOOOD+OO 

0.0000 

3 

-0. 1 OOOOOOOOD+O 1 

-0. 1 OOOOOOOOD+O 1 

0.0000 

4 

-0.300000000D+01 

-0 . 300000000D+0 1 

0.0000 

5 

0.30000000CD+01 

0 . 300000000D+0 1 

0.0000 


0. 15000000CD+01 

0. 1 50000000D+0 1 

0.0000 

7 

0. 10000000CD+01 

0. 1 OOOOOOOOD+O 1 

0.0000 

a 

0.249734643D+01 

0.225000000D+01 

11.0021 

9 

0.23000000CD+00 

0. 250000000D+00 

0.0000 

10 

-0. 17049343BB+00 

0. OOOOOOOOOD+OO 

-0. 1705 

• « 

-0. 184979870D+01 

-0. 1 50000000D+0 1 

23.3199 

12 

-0 . 200000000D+0 1 

-0 . 200000000D+0 1 

0.0000 

13 

0 . 600000000 D+0 1 

0. 600000000D+01 

0.0000 

14 

0 . 400000000D+01 

0 . 400000000D+0 1 

0.0000 

15 

0.31 3720643D+0 1 

0. 300000000D+01 

4.3735 

16 

0. 1 OOOOOOOOD+O l 

0. 1 OOOOOOOOD+O 1 

0.0000 

17 

0. 1 10000000D+02 

0. 110000000D+02 

0.0000 

13 

0 . 950000000 D+01 

0. 930000000D+01 

0.0000 

19 

0.900000000D+01 

0 . 900000000D+0 1 

0.0000 

20 

0. 1049921 16D+02 

0. 102500000D+02 

2.4313 

21 

0. 825000000D+01 

0 . 825000000D+0 1 

0.0000 

«■<» 

0.782284404 D+01 

0 . 800000000D+0 1 

-2.2144 

23 

0.615186643D+01 

0. 630000000D+01 

—3. 3559 

24 

0 . 600000000D+0 1 

0. 600000000D+01 

0.0000 


7. 

★ * * 
7. 

7. 

7. 

7. 

7. 

7. 

7. 

*** 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 


mcsh nr O'*- / a 
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** TEMPERATURE DISTRIBUTION ** MiTS* EL 


NODE NUMBER 


TEMPERATURES 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

i: 

12 

13 

14 

15 

16 


0 . 2 000000 OOOD+O 1 
O. 21 047978 18D-34 
1000000000D+01 
- . 30000 OOOOOD+O 1 
0. 30000 OOOOOD+O 1 
0. 10000 OOOOOD+O 1 
0. 7559057347D— 14 
2000000000D+0 1 
0. 600 00 OOOOOD+O 1 
0 . 4000000000D+0 1 
O. 3000000 OOOD+O 1 
0. 1000000000D+01 
0. 1 1 OOOOOOOOD+02 
0. 9000000000D+0 1 
0 . 80000000000+01 
0. 6000000000D+0 1 


( « /b/i) 
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TEMPERATURE DISTRIBUTION - MESH IV (8/8/8) 


NODE NUMBER 

NUMERICAL NUMBER 

ANALYTICAL VALUES 

ERROR 


1 

0 . 2000000000D+01 

0 . 200000000D+01 

0.0000 

% 

2 

0 . 2 1047978 18 D- 3 4 

0.000000000D+00 

0.0000 

*** 

3 

1000000000D+01 

-0. 100000000D+01 

0.0000 

% 

4 

3000000000D+01 

-0. 300000000D+01 

0.0000 

% 

5 

0. 3000000000D+01 

0 . 300000000D+01 

0.0000 

% 

6 

0. 1000000000D+01 

0. 100000000D+01 

0.0000 

% 

7 

0 . 7 5 59 0 57 54 7 D- 14 

0.000000000D+00 

0.0649 

*** 

8 

2000000000D+01 

-0. 200000000D+01 

0.0000 

% 

9 

0 . 6000000000 D+ 0 1 

0 • 600000000D+01 

0.0000 

% 

10 

0 . 4 000000000D+01 

0 . 400000000D+01 

0.0000 

% 

11 

0. 3000000000D+01 

0. 300000000D+01 

0.0000 

% 

12 

0. 1000000000D+01 

0. 100000000D+01 

0.0000 

% 

13 

0 . 1100000000D+02 

0. 110000000D+02 

0.0000 

% 

14 

0 . 9000000000D+01 

0 . 900000000D+01 

0.0000 

% 

15 

0. 8000000000D+01 

0 . 800000000D+01 

0.0000 

% 

16 

0 . 6000000000D+01 

0 . 6000000 00D+01 

0.0000 

% 


25 


ERROR ANALYSES 
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FINITE ELEMENT MESHES FOR ERROR ANALYSES I. 


Mesh I. 


Mesh II. 


Mesh III. 



1 12-noded element (max error - -0.1554) 

K 



8 12-noded elements (max error - -0.0475) 


64 12-noded elements (max error - -0.0119) 
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NODE 

NUMERICAL 

ANALYTICAL 

ERROR 

X ERROR 

1 

-0. 100000+01 

-0. 10000D+01 

0.0000 

0.0000 X 

2 

0. 23907D-33 

0. 000000+00 

0.0000 

**♦ 

3 

-0. 100000+01 

-0. 10000D+01 

0.0000 

0.0000 X 

4 

0.239070-33 

0. 000000+00 

0.0000 

*♦* 

3 

-0. 153440+00 

0. 000000+00 

-0.1334 

*•* 

6 

0.733160+00 

0. 730000+00 

-0.0168 

-2.2433 X 

7 

0. 10000D+01 

0. 100000+01 

0.0000 

0.0000 X 

8 

-O. 1 3344D+00 

0. OOOOOD+OO 

-0. 1334 

*** 

9 

0. 100000+01 

0. 10000D+01 

0.0000 

0.0000 X 

10 

-0. 133440+00 

0. 000000+00 

-0.1334 

*** 

11 

0 . 733 1 6D+00 

0. 73000D+00 

-0.0168 

-2.2453 X 

12 

0. 100000+01 

0. 10000D+01 

0.0000 

0.0000 X 



t (. t IZ n»df d 




** t\x error ■» -o.l55f- 
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Gj CO 10 U i.J U 10 Ut i.j t j f'j f'j f'j f'j fO f'.J IM fu I'J Pj 


* - ^ 

trICAL 

+N~i_'^TIC41. 

2 9 9 u 9 

X * 2.9 R On 

4 

10 GOOD +01 

-0. 1 OOCOD+Ol 

0. 0000 

0 . 0000 V. 


- 3 2 50C0D+CG 

-0 2 50C0D+00 

0. 0000 

0 . 0000 *. 


0. 20 5 23 £• -34 

0. OOCOOD+OO 

0. 0000 

-> 


- 0 1 C> GOOD "*“G l 

-0. 1 OOOOD+OI 

0. 0000 

0. 0000 7: 

f 

- 0 2 -5 000 D + GO 1 

-0. 2 5000D+00 

0. 0000 

0. 0000 7. 

c 

C 61057D-3+ 

0. OOOOOD+OO 

0. 0000 

■> ■> -> 

- 

-0 lOOCOD+Ol 

-0. 1 OOOOD+OI 

0 0000 

0. 0000 X 

; r* 

-0. 2 5 GOOD +00 

-0. 2 5000D+00 

0. 0000 

C. 0000 '/. 

9 

0 3C528D-34 

0. OOOOOD+OO 

0. 0000 

•*** 

1C 

-0 79747D+00 

-0. 75000D+00 

-0. 047 5 

-6. 3294 */. 

1 1 

-C 33249D+00 

-0. 3 1250D+Q0 

-0. 0200 

-6. 3931 */: 

4 \ 

i w 

-0 43250D-01 

0. OOOOOD+OO 

-0. 0432 

* 

4 ■— i 
1 - 

0. 18326D+00 

0. 1 8750D+00 

-0. 0042 

-2. 2618 7. 

14 

0. 2 5 GOOD +00 

0. 25000D+00 

0. 0000 

0. 0000 '/. 

1 5 

-0. 79747D+00 

-0. 75000D+00 

-0. 0475 

-6. 3294 7. 

1 6 

-0. 43250D-01 

0. OOOOOD+OO 

-0. 0432 


17 

0. 25000D+00 

0. 25000D+00 

0. 0000 

0. 0000 71 

18 

-0. 79747D+00 

-0. 75000D+00 

-0. 0475 

-6. 3294 7. 

19 

-0 33249D+00 

-0. 3 1 250D+00 

-0. 0200 

-6. 3981 7. 

20 

-0 43 2 SOD -01 

0. OOOOOD+OO 

-0. 0432 

■* 

2 1 

0 1S326D+0C 

0. 1 8750D+00 

-0. 004 2 

-2. 2618 7. 


0 2 5 GOOD +00 

0. 2 5000D+00 

0 0000 

0. 0000 71 

2 3 

-0 79 747D+00 

-0. 7 5000D+00 

-0. 047 5 

-6. 3294 7. 

2 4 

-0 43250D— 01 

0. OOOOOD+OO 

-0. 0432 

# *+ 

25 

C 25000D+00 

0. 2 5000D+00 

0. 0000 

0. 0000 7. 

2 b 

-0 79747D+00 

-0. 75000D+00 

-0. 0475 

-6. 3294 7. 

^ T 

zL / 

-0 33249D+00 

-0. 3 1250D+00 

-0. 0200 

-6. 3981 7. 

-v 0 

& w 

-0 43250D— 01 

0. OOOOOD+OO 

-0. 0432 


29 

0 18326D+00 

0. 1 8750D+00 

-0. 0042 

-2. 2618 7. 

2 0 

Q. 2 r GOOD +00 

0. 2 5000D+00 

0. 0000 

0. 0000 7, 

7 - 

0 1 57470-0: 

0. OOOOOD+OO 

0. 0157 


b — 

0. 77 1 23D+00 

0. 7 50C0D+00 

0. 0212 

2. 8301 71 

T 1 ‘T 

0. 1 OOOOD+Oi 

0 1 OOOOD+OI 

0. 0000 

0. 0000 7. 

34 

0 1 5747D— 01 

0. OOOOOD+OO 

0 0157 

■Jf-if- *> 

3 5 

C 771 23D+00 

0. 7 5000D+00 

0. 0212 

2. 9301 X 

3 b 

0 10 GOOD +01 

0. 1 OCOOD+Ol 

0 0000 

0. 0000 X 

2 7 

0 1 3 747D-01 

0. OOOOOD+OO 

0. 0157 


39 

0. 771 23D+00 

0. 75000D+00 

0 0212 

2. 8301 X 

3 9 

0 1 OOOOD+OI 

G 1 OOOOD+OI 

0 0000 

0. 0000 X 


MCI 14 IT % 

12 

C ffror 

* - 0. 0^.75) 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


r»i ro 


ANAL Y f 


A i 

r-._ 


10 
1 1 


1 7 
1 0 
1 ? 





-0 250000+00 
— 0 “• 2 5 00 D — 0 L 
0 '’43210-35 
25 GOOD +00 
-C *250CD-0i 
0. 381610-35 
-0 19937D+00 
-0. 331230-01 
-0 108120-01 
... 458150-01 
0. 625000-01 
-0 199370+00 
-0 1081 2 D -01 
0. 62500D-01 
-0. 199370+00 
-0 83123D-01 
-0 10812D-01 
0 45S15D-C1 
C 425003*— Ol 
-0. 1 99370+00 
-0 1O012D-O1 
0 425000-01 
-0. 199370+00 
-C S3 1230-01 

-0 10ei2D-01 

0. 458150-01 
0 625000-01 
0. 39 3680-02 
0. 192810+00 
0 250000+00 
C 3°36SD-02 
0 1 0 28 ID +00 
0. 250000+00 
C 393680-02 
0. 192810+00 
0 250000+00 


-0. 2 5000D+00 
-0 6 2 SCO D -01 
0. 000000+00 
-0. 25000D+00 
-0. 625000—01 
0. OOOCOD+OO 
-0. 250000+00 
-0. 62500 D-01 
0. OOOOOD+GO 
-0. 1 8750D+00 
-0. 78125D-01 
0. OOOOOD+OO 
0. 46875 D-01 
0. 625000-01 
-0 1 8750D+00 
0. OOOOOD+OO 
0. 62500D-01 
-0. 187500+00 
-0. 781250-01 
0. OOOOOD+OO 
0. 468750-01 
0. 62500D-01 
-0. 187500+00 
0. OOOOOD+OO 
0. 62500D-01 
-0. 187500+00 
-0. 781250-01 
0. OOOOOD+OO 
0. 46875D-01 
0. 625C0D-01 
0. OOOOOD+OO 
0. 1 8750D+00 
0. 25000D+00 
0. OOOOOD+OO 
0 187500+00 
0. 250C0D+00 
0. OOOOOD+OO 
0. 187500+00 
0. 25000D+00 



0 . 0000 
0 . 0000 
0 . 0000 
0 . 0000 
0 . 0000 
0 . 0000 
0 . 0000 
0 . 0000 
0 . 0000 
- 0 . 0119 
- 0. 0050 
- 0 . 0103 
- 0. 001 1 
0 . 0000 
- 0 . 0119 
- 0 . 0108 
0 . 0000 
- 0. 0119 
- 0. 0050 
- 0. 0108 
- 0 . 001 1 
0 . 0000 
- 0 . 0119 
- 0. 0108 
0 . 0000 
- 0. 0119 
- 0 . 0050 
-0 0108 
- 0 . 001 1 
0 0000 
0 . 0039 
0 . 0053 
0 . 0000 
0 . 0039 
0 . 0053 
0 . 0000 
0 . 0039 
0 . 0053 
0 . 0000 


\ ERRCF 


0. 0000 
0. 0000 
-<*■->-> 

0. 0000 
0. 0000 

0. 0000 

C . 0000 
•*** 

- 6 . 3294 
- 4 . 3981 
** -> 

- 2 . 2613 
0 . 0000 
- 6 . 3294 

0 . 0000 
- 6 . 3294 
- 6 . 3981 
*** 

- 2 . 2618 
0 . 0000 
- 6 . 3294 

0 . 0000 
- 6 . 3294 
- 6 . 3981 
**» 

- 2 . 2618 

0 . 0000 

2 . 8301 
0 . 0000 

2 . 8301 
0 . 0000 
*■*-* 

2 . 8301 
0 . 0000 


V. 


■/ 

/• 

•/ 



•/ 

»/ 

/• 


7 . 

0/ 

A 

7. 


7. 


7. 


a/ 

A 

7 . 
• ■ 
/« 


7. 


7. 


7. 

7. 


V 

V 


V 

7 . 
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ERROR ANALYSES II. (1 element, 8 element solutions) 


Boundary Conditions 

ti *- 0 - 3 y 2 + 4zZ 

T|y_ 0 - 2x 2 + 4z 2 
Tl^o - 2x 2 + 3y 2 



fill - 4 

ax i^i 

fill - 6 

dyly«l 

fill - 8 

3z | Z-1 


I.® 



V 2 T-18 
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FINITE ELEMENT MESHES FOR ERROR ANALYSES II. 



Mesh I. 1 12-noded element (max error — 0.983) 



Mesh II. 8 12-noded elements (max error — 0.2733) 


32 




O«000Njthoi4*UW 


NODE NUMERICAL 


1 0. 20000D+01 

-0. 20692D-32 
0. 60000D+01 
0. 40000D+01 
0. 50000D+01 
0. 35000D+01 
0. 30000D+01 
0. 61 973D+01 
0. 40000D+01 
0. 99830D+01 

11 0. 79362D+01 

12 0. 70000D+01 


ANALYTICAL 


0. 20000D+01 
0. OOOOOD+OO 
0. 60000 D+01 
0. 40000 D+01 
0. 50000D+01 
0. 35000D+01 
0. 30000D+01 
0. 60000D+01 
0. 40000D+01 
0. 90000D+01 
0. 75000D+01 
0. 70000D+01 


CrJC-ort. 3 

MctH X C 


ERROR % ERROR 


0. 0000 

0. 0000 

7. 

0. 0000 

*** 


0.0000 

0. 0000 

7. 

0.0000 

0. 0000 

7. 

0. 0000 

0. 0000 

7. 

0. 0000 

0. 0000 

7. 

0. 0000 

0. 0000 

7. 

0. 1973 

3. 2880 

7. 

0. 0000 

0. 0000 

7. 

0. 9830 

10. 9223 

7. 

0. 4362 

5. 8162 

7. 

0. 0000 

0. 0000 

7. 


e^rroy" - ° 


33 


NODE 

NUMERICAL 

ANALYTICAL 

ERROR 

7. ERROR 

1 

0. 200000+01 

0.200000+01 

0.0000 

0.0000 X 

2 

0.50000D+00 

0. 500000+00 

0.0000 

0.0000 % 

3 

-0. 245660-33 

0. 00000D+00 

0.0000 

*♦* 

4 

0.300000+01 

0.300000+01 

0.0000 

0.0000 X 

5 

0. 15000D+01 

0. 15000D+01 

0.0000 

0.0000 X 

6 

0. 10000D+01 

0. 10000D+01 

0.0000 

0.0000 7. 

7 

0.600000+01 

0.600000+01 

0.0000 

0.0000 7. 

8 

0.450000+01 

0. 45000D+01 

0.0000 

0.0000 X 

9 

0. 400COD+01 

0. 40000D+01 

0.0000 

0.0000 7. 

10 

0.27500D+01 

0. 27500D+O1 

0.0000 

0.0000 7. 

11 

0. 187500+01 

0. 18750D+01 

0.0000 

0.0000 X 

12 

0. 125G0D+01 

0. 12500D+P1 

0.0000 

0.0000 % 

13 

0.37300D+00 

0. 87500D+00 

0.0000 

0.0000 X 

14 

0. 75C00D+00 

0.750000+00 

0.0000 

0.0000 X 

15 

0.30766D+01 

0. 30000D+01 

0. 0766 

2.5547 X 

16 

0. 156320+01 

0. 15000D+01 

0. 0632 

4.2136 7. 

17 

0. 1 00000+01 

0. 100000+01 

0. 0000 

0.0000 7. 

18 

0.400720+01 

0. 37500D+01 

0. 2572 

6.8592 7. 

19 

0. 307490+01 

0. 2B750D+01 

0. 1999 

6.9525 7. 

20 

0. 251B3D+01 

0. 22500D+01 

0. 2683 

11. 9239 X 

21 

0. 19737D+01 

0. 18750D+01 

0. 1187 

6.3326 X 

22 

0. 17500D+01 

0. 17500D+01 

0. 0000 

0.0000 X 

23 

0.51 63SD+0 1 

0.500000+01 

0. 1688 

3.3758 7. 

24 

0.36437D+01 

0. 35000D+01 

0. 1437 

4.1055 7. 

25 

0.300000+01 

0. 30000D+01 

0. 0000 

0.0000 7. 

26 

0. 701 120+01 

0. 67500D+01 

0.2612 

3.8692 X 

27 

0.60*360+01 

0. 58750D+01 

0.2236 

3.8064 7. 

28 

0. 552330+01 

0. 52500D+01 

0. 2733 

5.2049 X 

29 

0. 500420+01 

0.487500+01 

0. 1292 

2.6502 X 

30 

0. 475000+01 

0.475000+01 

0. 0000 

0.0000 7. 

31 

0.500000+01 

0. 50000D+01 

0.0000 

0.0000 X 

32 

0.35000D+01 

0.350000+01 

0.0000 

0.0000 X 

33 

0.300000+01 

0. 30000D+01 

0.0000 

0.0000 7. 

34 

0.5S631D+01 

0. 60000D+01 

-0.1319 

-2.1991 7. 

35 

0. 43490D+01 

0. 45000D+01 

-0. 1510 

-3.3546 7. 

36 

O.4.>00OD+Ol 

0. 40000D+01 

0.0000 

0.0000 7. 

37 

0.333380+01 

0. 90000D+01 

-0.1112 

-1.2355 7. 

38 

0. 73743D+01 

0. 75000D+01 

-0. 1257 

-1.6756 7. 

39 

0.7CC00D+01 

0. 70000D+01 

0. 0000 

0.0000 7. 


tttftort- AMrxt'r \ n IT 

I ( c»"r*v 

\2 noi«J ele«\o%4-s 
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V.b. Structural Analysis Applications 
DISPLACEMENT ANALYSIS (END LOADING) 



Boundary Conditions 

at x-0 built-in end 

at all points w is prescribed 

The problem is solved by using different finite element mesh types. 
Results obtained demonstrate the power of the mesh interface elements. 
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FINITE ELEMENT MESHES FOR END LOADING PROBLEM 


ll 


/T . A S. T? 

k - . - ‘ 

* / 

( 

r ? 

- ^ * * 

f ^ 

m 

• 

i ^ 

r 

y- i 

a 


Mesh I. 3 12-noded elements 



Mesh II. 3 20-noded elements 


Ifc 


z 

L 7 

' ^7' 

i 

1 

• —4 
* 

* — 


BBS 


Mesh III. 3 8-noded elements 



Mesh IV. 8 -noded/12 -noded/20-noded elements 
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TEST CA5E:SbAM 3 »/li/*o -ELEM < ****** 

10 


26 

26 

3 

3 3 

28 

1 

4 

8 6 

1 

3 

0 0 

8 











1 

1 

14 

1 

3 

6 

4 

2 

6 

7 

3 

12 
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1 

13 

3 

9 

13 

8 

6 

11 

13 

7 




16 

10 

14 

12 





20 











3 

1 

1 6 

9 

21 

27 

13 

11 

23 

23 

13 




17 

ia 

20 

16 

10 

22 

26 

14 




18 

24 

19 

12 






1 




0.0 


0.0 

2 




0.0 


1.0 

3 


r - m 


1.0 


1.0 

4 


.* , 


1.0 


0.0 

5 


1 . 


0.0 


0.0 

6 


1 . 


0.0 


1.0 

7 


1 . 


1.0 


1.0 

8 


1 . 


1.0 


0.0 

9 


4 , 


0.0 


0.0 

10 


*' > 


0.0 


0.5 

1 1 


2 . 


0.0 


1.0 

12 


4 . 


0.3 


1.0 

13 


* • 


1.0 


1.0 

14 


*’ > 


1.0 


0.5 

15 


* . 


1.0 


0.0 

16 


2‘ 0 


0.3 


0.0 

17 


3 . 


0.0 


0.0 

18 


3 . 


0.0 


1.0 

19 


3 . 


1.0 


1.0 

20 


3 . 


1.0 


0.0 * 

21 


3 . 


0.0 


0.0 

22 


3 . 


0.0 


0.5 

23 


3 . 


0.0 


1.0 

24 


3 . 


0.3 


1.0 

23 


3 . 


1.0 


1.0 

26 


3 . 


1.0 


0.5 

27 


3 . 


1.0 


0.0 

28 


3 , 


0.3 


0.0 

1 

2 . 

9 *E +07 

0.3 


. 0.0 

1 

1 


1 

0 . 

1 

0 . 

2 

1 

*• # 

1 

0 . 

1 

0 . 

3 

1 

f - # 

1 

0 . 

1 

0 . 

4 

1 

*• # 

1 

0 . 

1 

0 . 

3 

0 

r m 

c* 

0 . 

1 

0 . 

6 

0 

*: m 

o 

0 . 

1 

o. 

7 

0 

0 

0 

0 . 

1 

0 . 

0 

0 


0 

0 . 

1 

0 . 

9 

0 


0 

0 . 

1 

0 . 

10 

o 

m 

c* 

o. 

1 

0 . 

11 

0 


c» 

0 . 

1 

o. 

12 

0 

'• # 

0 

0 . 

1 

o. 

13 

0 


0 

0 . 

1 

0 . 

14 

0 


0 

0 . 

1 

o. 

13 

0 

.'i ^ 

0 

0 . 

1 

0 . 

16 

0 


0 

0 . 

1 

0 . 

17 

0 

.*« 

0 

0 . 

1 

0 . 

18 

0 

r , # 

o 

0 . 

1 

0 . 

19 

0 


o 

0 . 

1 

0 . 

20 

0 


0 

0 . 

1 

0 . 

21 

0 

0 « 

0 

0 . 

1 

0 . 

22 

0 

v. 

0 

0 . 

1 

0 . 

23 

0 

0 , 

0 

0 . 

1 

0 . 

24 

0 

0 * 

0 

0 . 

1 

o. 

23 

0 

V . 

0 

0 . 

1 

0 . 

26 

0 

V « 

0 

o. 

1 

o. 

27 

0 

V % 

0 

0 . 

1 

o. 

28 

0 

0 , 

0 

o. 

1 

0 . 

21 


0 « 


30 v. 

0 . 


22 




300 . 

0 . 


23 


v . 


300 . 

0 . 


24 


0 , 


300 . 

o. 


23 


s 


300 . 

0 . 


26 




so;-. 

0 . 


27 


0* 


so:. 

0 . 


28 


V * 


so:. 

0 . 



37 


DEFLECTIONS AT THE TIP (END LOADING) 
(*“12 , y-1, z-0) 


Mesh Type 

u 

V 

w 

Mesh I 

-0.002388 

0.010315 

Presc . 

Mesh II 

-0.00323 

0.01373 

Presc . 

Mesh III 

-0.00211 

0.00928 

Presc . 

Mesh IV 

-0.00236 

0.00983 

Presc. 
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DISPLACEMENT ANALYSIS (AXIAL LOAD) 
Mesh I (8 12-noded elements) 



Node Tip Deflections 

Number u y 


19 

20 
21 
22 

23 

24 

25 

26 
27 


0.2347 x 10' 3 
0.1242 x 10' 3 
0.2347 x 10' 3 
0.1242 x 10‘ 3 
0.9639 x 10* 4 
0.12417 x 10' 3 
0.2347 x 10' 3 
0.1242 x 10' 3 
0.2347 x 10' 3 


0.4313 x 10’ 4 
0.2760 x 10 -4 
0.4313 x 10* 4 
0.0 
0.0 
0.0 

*0.4313 x 10’ 4 
-0.2760 x 10' 4 
-0.4313 x 10’ 4 


0.4313 x 10' 4 

0.0 

-0.4313 x 10' 4 
0.2760 x 10' 4 
0.0 

-0.2760 x 10' 4 
0.4313 x 10' 4 
0.0 

-0.4313 x 10* 4 
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DISPLACEMENT ANALYSIS (AXIAL LOAD) 
Mesh II (8 20-noded elements) 



Node 




T1d Deflections 



u 


V 

w 

19 

0.5233 

X 

10' 3 

0.1393 x 10' 3 

0.1393 x 10" 3 

20 

0.2191 

X 

10' 3 

0.505 x 10"* 

0.0 

21 

0.5233 

X 

10' 3 

0.1393 x 10" 3 

-0.1393 x 10" 3 

22 

0.2191 

X 

10' 3 

0.0 

0.505 x 10'* 

23 

0.9122 

X 

10'* 

0.0 

0.0 

24 

0.2191 

X 

10' 3 

0.0 

-0.505 x 10'* 

25 

0.5233 

X 

10' 3 

-0.1393 x 10" 3 

0.1393 x 10' 3 

26 

0.2191 

X 

10' 3 

-0.505 x 10"* 

0.0 

27 

0.5233 

X 

10* 3 

-0.1393 x 10" 3 

-0.1393 x 10' 3 
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DISPLACEMENT ANALYSIS (AXIAL LOAD) 
Mesh III (8 8- node d elements) 


n 



Node Tip Deflections 

Number u y 


19 

0.2265 

X 

10' 3 

20 

0.1252 

X 

10' 3 

21 

0.2265 

X 

10' 3 

22 

0.1252 

X 

10' 3 

23 

0.9994 

X 

10'* 

24 

0.1252 

X 

10' 3 

25 

0.2265 

X 

10' 3 

26 

0.1252 

X 

10' 3 

27 

0.2265 

X 

10' 3 


0.4194 x 10'* 
0.2727 x 10'* 
0.4194 x 10'* 
0.0 
0.0 
0.0 

-0.4194 x 10'* 
-0.2727 x 10"* 
-0.4194 x 10'* 


0.4194 x 10'* 

0.0 

-0.4194 x 10'* 
0.2727 x 10'* 
0.0 

-0.2727 x 10'* 
0.4194 x 10'* 
0.0 

-0.4194 x 10'* 
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DISPLACEMENT ANALYSIS (AXIAL LOAD) 
Mesh IV (4 8-noded/4 12-noded elements) 



Node 




TId Deflections 



u 


V 

V 

19 

0.4129 

X 

10' 3 

0.722 x 10'* 

0.722 x 10'* 

20 

0.2142 

X 

10' 3 

0.3731 x 10'* 

0.0 

21 

0.4129 

X 

10‘ 3 

0.722 x 10"* 

-0.722 x 10'* 

22 

0.2142 

X 

10' 3 

0.0 

0.3731 x 10‘* 

23 

0.1089 

X 

10" 3 

0.0 

0.0 

24 

0.2142 

X 

10' 3 

0.0 

-0.3731 x 10"* 

25 

0.4129 

X 

10' 3 

-0.722 x 10'* 

0.722 x 10‘* 

26 

0.2141 

X 

10' 3 

-0.3731 x 10'* 

0.0 

27 

0.4129 

X 

10' 3 

-0.7221 x 10'* 

-0.722 x 10'* 
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APPENDIX 

C 

C 

C 

IMPLICIT REAL*8(A-H, 0-Z ) 

REAL*8 LU 
C 

DIMENSION COOR (200* 3) * NCONEC (20, 21 ),NB0UN(100, 8), NBKIND(100, 3) 
DIMENSION CBOUN< 100, 3) , GST IF (200, 200) , GL0AD(200), GTEMP <200 ) 
DIMENSION STIF (20, 20), EL0AD(8) , ESTIF(8, 8). TEMP (200 ) 

DIMENSION INDEXR(200), LU(200, 200) . W(6) , Z(6 ), SHAP(200) 

DIMENSION TEMP AN(200), C (200), FORCE (20), BL0AD(20) 

C 

C N1 = Max number of nodes 

C N2= Max number of boundary surfaces 

C N3= Max number of nodes in any element 'in this case N3=20 ' 

C N4= Max number of elements 

C 

DATA N1,N2,N3, N4/200, 100,20,20/ 

C DATA ( W(M ) , M=l, 3)/. 88888B88B888B9. . 55553555555356, . 35555555555556/ 

C DATA ( Z ( M ) , M= 1 , 3 ) / . 0, . 774596669241483, -.774596669241483/ 

C 

DATA (W(M),M=1,6>/. 467913934572691, .467913934572691, 

> . 360761573048139, . 360761573048139, 

> . 171324492379170, . 171324492379170/ 

C 

DATA ( Z (M) , M=1 , 6)/. 238619186083197,-. 238619186083197, 

> . 661209386466265, -. 661209386466265, 

> . 932469514203152, -. 9324695142031 52/ 

C 

DATA PENALT/ 100000000000000000000000000000000. / 

C 

C 

OPEN (UNIT=1,FILE='0NUR. DAT STATUS='OLD', READOM.Y ) 

OPEN (UNIT=2, FILE= / ONUR. OUT ', STATUS* 'NEW' ) 

C 

C 

c 

C READ NUMBER OF NODES, ELEMENTS, BOUNDARIES 

C 

READ(1,*) NNODES 
READ<1,*> NELEM 
READd,*) NBOUND 
C 

C READ COORDINATES 

C 

C 

DO 10 1 = 1, NNODES 

READd, *) COOR(I, 1 ) , COOR (I, 2),C00R(I,3) 

10 CONTINUE 
C 
C 

C READ CONNECTIVITIES 

C 

DO 20 1=1, NELEM 
C 

C * NCONEC ( I , 1 ) defines the type of the element '8,12 or 20' 

C 

READd, 21) NCONEC < 1 , 1) 

READd, 21) (NCONEC < I, J>, U=2, NCONEC (I. 1)+1) 

READd, *) FORCE ( I ) 
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20 CONTINUE 

21 FORMAT <2013) 


READ BOUNDARY PARAMETERS 

Surface number* node number* element number 
DO 25 1=1* NBOUND 

READ<1.*> NBKINDd, 1 ), NBKINDd* 2>» NBKINDd* 3) 
25 CONTINUE 


READ BOUNDARY CONNECTIVITIES 


DO 30 1=1, NBOUND 

READ< 1* 31)' <NBOUN( I , J) * J=l* NBKINDd. 2)) 

30 CONTINUE 

31 FORMAT < 81 3 ) 


READ BOUNDARY CONSTANTS 

Conduction coef. * convection coef. 

DO 40 1 = 1* NBOUND 

READ (1**) CBOUN d * 1 ) , C BOUN (I * 2 ) , CBOUN (1,3) 
40 CONTINUE 


READ THE CONDUCTION COEFFICIENT 


OF THE DOMAIN 


READd,*) COEF 

PRESPECIFIED NODAL TEMPERATURES 

DO 43 1=1, N1 
TEMP < I )=0. DO 
43 CONTINUE 


* NTEMP= Number of nodes where temperature is aiready specified 

* ITEMP= Node number 


DO 44 1=1* N 1 

GLOAD( I )=0. DO 
CTEMP ( I )=0. DO 
DO 44 11=1, N1 

GST I F ( I * II >=0. DO 
44 CONTINUE 


READd,*) NTEMP 
DO 45 1 = 1, NTEMP 
READ(1,*> ITEMP, TEMP( I TEMP > 
GSTIF ( ITEMP, ITEMP ) =PENALT 
45 CONTINUE 
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c 

c 

c 

READ<1.*> Cl. C2. C3. C4. C5, C6.C7 

ri i i i i i i i i i i i • i t i i i i i 1 1 i i i i i i 1 1 i i 1 1 i i i .) i i 1 1 • i i ■ i i i • t i 

C 

DO 46 1=1 . NNODES 

7EMPAN ( I)=C1+ COOR ( I. 1 ) *<C2+C5*C00R ( I. 1)> 

> + COOR < 1.2) *<C3+C6*C00R( 1.2) ) 

> + COOR < 1.3) *<C4+C7#C00R< 1,3) ) 

46 CONTINUE 

£i * i « i • ■ * i i i i i i i i { i i i i [• i 1 1 i i j i j i 1 1 i i 1 1 i { 1 1 i i f ; i i 1 1 i ■ m [ 

C 

DO 50 1=1. N1 

CLC AD < I ) =PENALT *TEMP ( I ) 

OTEMP < I )=0. DO 
50 CONTINUE 
C 

c 

c 

c 

C CALCULATE THE GLOBAL STIFFNESS MATRIX 
C 
C 
C 

DO 100 1=1 . NELEM 
C 

CALL ELEM ( COOR . NCONEC . STIF. COEF. I, Nl. NS. N3. N4. W. Z, BLOAD. 

> FORCE. SHAP) 

C 

C 

C 

C 

C 

DO 60 J=l. NCONEC ( 1 . 1) 

IY=NCONEC< I, J+l) 

GLOAD( IY)=GLOAD( IY )+BLOAD( J) 

C 

DO 59 K=1 NCONEC < 1 . 1) 

C 

IX=NCONEC< I.K+1) 

C 

GSTIF < IY. IX > = GST IF < IY. IX > + STIF ( J, K / 

C 

c 

59 CONTINUE 

60 CONTINUE 
100 CONTINUE 

C 

C CALCULATE THE GLOBAL LOAD VECTOR 

C 

DO 110 1=1 . NBOUND 

N=NBKIND( I. 1) 

INODE=NBKIND< I. 2) 

IELEM=NBKIND( I, 3) 

C 

C 

c 

IF (N. EQ. 1) CALL BOUNDR ( 0. DO. COOR. NBOUN. NCONEC. ELOAD, IELEM, INODE. I. 

> W, Z.N1.N2. N3.N4. SHAP.ESTIF) 
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IF (N. EQ. 2) 

> 

IF (N. EQ. 3) 

> 

IF (N. EQ. 4) 

> 

IF (N. EQ. 5) 

> 

IF (N. EQ. 6) 

> 


CALL BOUNDR ( 1. DO, COOR* NBOUN, NCONEC, ELQAD, IELEM, INODE, I, 
U, Z, Nl, N2, N3. N4» SHAP, ESTIF) 

CALL BOUNDS (0. D0,C00R. NBOUN, NCONEC, ELOAD. IELEM, INODE, I, 
W, Z, Nl, N2, N3, N4, SHAP, ESTIF) 

CALL BOUNDS (1. DO, COOR, NBOUN, NCONEC, ELOAD, IELEM, INODE, I, 
W, Z, Nl, N2. N3.N4, SHAP. ESTIF) 

CALL BOUNDT (0. DO, COOR, NBOUN. NCONEC, ELOAD, IELEM, INODE, I, 
W, Z, Nl, N2, N3. N4» SHAP, ESTIF) 

CALL BOUNDT (1. DO, COOR, NBOUN, NCONEC, ELQAD, IELEM, INODE, I, 
U, Z, Nl, N2, N3.N4. SHAP. ESTIF) 


DO 105 J=l, INODE 
LL=NBOUN(I, J) 

GLOAD(LL)=GLOAD(LL)+<CBOUN(I. 1 )+CB0UN( I, 2)*CB0UN< 1,3) )*EL0AD< J) 
DO 105 K=l, INODE 
LK=NBOUN(I, K) 

OSTIF (LL, LK)=GSTIF(LL» LK) +CBOUN( I. 2)*ESTIF< J, K) 

105 CONTINUE 
110 CONTINUE 


CALL DECOMP (NNODES, GSTIF, LU, INDEXR.Nl) 

CALL SOLVE (NNODES, LU, GLOAD, GTEMP. INDEXR, Nl ) 


i i i i i ft i i i i i i i i i i i i * i ii i ! ! \ 


i i i 1 1 i i i i i i i i i i i i i i i i * • • i j • ; j M!!!!! ! 


WR ITE ( 2, 999) 

DO 99 S K=l. NNODES 
C (K>=( GTEMP <K) -TEMPAN( K) ) 

IF (DABS (TEMP AN (K) ) . LT. 0.00000000000001) THEN 
WRITE (2, 1000) K, GTEMP (K), TEMPAN(K). C (K) 

ELSE 

CK=C (K)*100/DABS (TEMPAN(K) ) 

WRITE(2, 1001) K, GTEMP (K) , TEMPAN(K). C (K), CK 
ENDIF 

998 CONTINUE 

999 FORMAT ( '1 ',//, 2X, ' NODE 5X. ' NUMERICAL 5X, 

>' ANALYTICAL ' , 5X. * ERROR '» 5X» * X ERROR '» /» 2X, 

> ' ', 5X, ' ', 5X, 

>, • , 5X, ' ', 5X, ' ’.//) 

1000 FORMAT (3X , 13, 6X. D12. 5, 5X, D12. 5, 4X. F9. 4. 8X. * ***') 

1001 FORMAT (3X, 13, 6X, D12. 5, 5X, D12. 5.4X.F9. 4, 5X.F9. 4. ' V.') 


STOP 

END 
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C %%%%%%%X%%%XXXXXXXXXXXXXXXXXZXXXXXXXXXXXXX%XXXXXXXX%XXXXXXXXXXXXXXXX%XX 
c 


SUBROUTINE ELEM(COOR, NCONEC, ST IF# COEF* IJK, Nl. N2* N3« N4> W, Z, 
> BLOAD, FORCE, SHAP) 

C 

C THIS PROGRAM CALCULATES THE STIFFNES MATRIX FOR 8/12/20 
C NODED FINITE ELEMENTS 
C 


IMPLICIT REAL*8(A-H» 0-2) 

DIMENSION COOR (Nl» 3), W(6), Z (6) , STIF<N3. N3) . SHAPE(3, 20) 
DIMENSION SJ(3. 3), DJ(3,3), ST(3, 20) , NCONEC ( N4, N3+1 ) 
DIMENSION BLOAD (20), FORCE ( N4), SHAP(Nl) 

C ####### ###################### 

DIMENSION TRAN (3, 3) 

C### ######## #-t################ 

c 


R=0. DO 
S=0. DO 
T=0 DO 


DO 20 1 = 1 , NCONEC ( I JK, 1 ) 

BLOAD (I )=0. DO 
DO 20 J=l, NCONEC (IJK, 1 > 

STIF < I, J)=0. DO 
20 CONTINUE 
C 

C ############################################## 
CALL TRANS ( IJK. Nl, N3. N4, COOR, NCONEC, TRAN) 
C### #################################### ####### 


DO 1000 1 — 1,6 

R=0. 5D0*( 1 . DO+Z ( I > ) 

DO 900 J=l. 6 

S=0. 5D0*< 1. DO+Z ( J> ) 

DO 800 K= 1 , 6 

T=0. 5D0*( 1. DO+Z (K) ) 

IF (NCONEC ( J JK, 1 ). EQ. 8) CALL ELEM08(R, S. T, SHAPE) 
IF (NCONEC ( IJK, 1). EQ. 12) CALL ELEM12(R, S, T, SHAPE) 
IF (NCONEC < IJK, 1>.EQ. 20) CALL ELEM20(R, S, T, SHAPE) 


IF(FORCE( I JK). EQ. 0. DO) GOTO 29 
NC =NCONEC < IJK, 1 ) 

IF(NC. EQ. 8> CALL SHAP08(R, S, T, NCONEC. SHAP. Nl. N3, N4, IJK) 
IF(NC.EQ. 12) CALL SHAP 12(R, S. T. NCONEC, SHAP. Nl. N3, N4, IJK) 
IF(NC. EQ. 20) CALL SHAP20(R, S. T. NCONEC, SHAP. Nl. N3. N4, IJK) 


29 DO 30 11=1, 3 
DO 30 JJ=1 , 3 

SJ(II, JJ)=0. DO 
DO 30 KK=1 , NCONEC ( IJK, 1 > 

IA=NCONEC ( IJK, KK+1 ) 

C 

SJ( I I, JJ)=SJ( II, J J ) + SHAP E ( 1 1 . KK ) *CGOR ( I A , Jj> 

30 CONTINUE 
C 
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c ######################################### 

CALL MULTI P (TRAN. SJ) 

C ##############*########################## 

C 

DET=SJ( 1. 1 )*(SJ(2, 2)*SJ(3» 3>-SJ(3. 2)*SJ(2. 3) ) 
DET=DET+SJ ( 1. 3 )*<SJ(2. 1 )*SJ(3, 2)-SJ(3» 1 )*SJ(2, 2) ) 
DET=DET-SJ( 1' 2)*<SJ(2» 1 )*SJ(3, 3)-SJ<3, 1 )*SJ(2, 3) ) 


DJ< 1. 1 ) = (SJ(2. 2)*SJ(3, 3)-SJ(3, 2)*SJ<2. 3>>/DET 
DJ (2. 2 )=(SJ( 1. 1)*SJ(3. 3)-SJ<3, 1)*SJ<1» 3>>/DET 
DJ(3. 3 ) = <SJ( 1» 1 )*SJ(2, 2)-SJ< 1, 2)*SJ<2» 1))/DET 
DJ< 1. 2)=(SJ(3» 2)*SJ(1, 3)-SJ( 1. 2)*SJ<3» 3))/DET 
DJ< 1* 3> = <SJ(i# 2)*SJ<2» 3)-SJ<2, 2)*SJ(1» 3) )/DET 
DJ (2* 1 )=(SJ(3» 1 )*SJ(2, 3)-SJ(2. 1>*SJ<3. 3>>/DET 
DJ(2. 3> = (SJ(2. 1)*SJ(1, 3)-SJ(l* i)*SJ(2* 3>)/DET 
DJ (3. 1 )=(SJ(2» 1 )*SJ(3» 2)-SJ(3» i)*SJ<2. 2))/DET 
DJ (3. 2)=(SJ(3. 1)*SJ(1, 2)-SJ< 1» 1 )*SJ(3» 2>>/DET 


DO 40 111 = 1,3 

DO 40 JJJ= 1, NCONEC ( IJK, 1 > 

ST(I II, JJJ>=0. DO 
DO 40 KKK= 1 , 3 

ST( III, JJJ)=ST( III, JJJ)+DJ( III# KKK)*SHAP£(KKK, JJJ) 
40 CONTINUE 


DO 50 11=1, NCONEC ( IJK, 1) 

IRZ=NCONEC ( IJK, 11 + 1) 

BLOAD< 1 1 )=BL0AD(I1 >+W( I )*W ( J)*W(K)*FORCE( I JK)*£HAP ( IRZ )* 
> DET/8. DO 

DO 50 Jl = l, NCONEC ( IJK, 1) 

ALPHA=0. DO 
DO 60 Kl=l , 3 

ALPHA=ALPHA+ST ( K 1 , II >*ST(Kl. Jl> 

60 CONTINUE 

ALPHA=ALPHA*COEF*DET/B. DO 

STIF (II, J1 )=STIF (II, J1)+W<I)*W(J)*W(K)*AlPHA 

50 CONTINUE 
800 CONTINUE 
900 CONTINUE 
1000 CONTINUE 


RETURN 

END 


C XXXXXXXXX%XX7ai%XXX7lXXXXX%XXXXKXX%XXXXXXXXXXXXXXXXX%X%XXXXXXXXX > X 


SUBROUTINE BOUNDR ( R, COOR. NBOUN, NCONEC, ELOAD, IJK, INODE, IB, 

W, Z, Nl, N2, N3. N4. SHAP, ESTIF) 


IMPLICIT REAL+8 ( A-H, 0-Z ) 

DIMENSION COOR (Nl. 3), W(6>. Z (6 ) . SHAPE(3. 20) 
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DIMENSION S J ( 3 , 3 ) , SHAP (Nl ), NB0UN(N2» 8) , NCONEC <N4, N3+1 ) 
DIMENSION ELOAD (8) » ESTIF(B, 8) 

C ################### ############### 

DIMENSION TRAN (3, 3) 

C ft#######*#######**#############*## 

C 

Ml -NCONEC ( IJK, 1) 

DO 10 1-1,8 
ELOAD ( I )=0. DO 
DO 10 I 1-1 i 8 
ESTIF(I, II >=0. DO 
10 CONTINUE 
C 

C #####################*##*#####«###«#*######### 

CALL TRANS ( IUK# Nl. N3. N4. COOR, NCONEC. TRAN) 
C################### ##################»#«#«#### 

c 

c 

DO 50 J= 1,6 

S=0. 5D0*( 1 . D0+2 ( J) > 

DO 40 K= 1,6 

T=0. 5D0*( 1 . DO+Z (K) ) 

C 

C 

IF (Ml. EQ. 8) CALL SHAP08 (R# S, T, NCONEC, SHAP, Nl, N3, N4. IJK) 
IF (Ml. EQ. 8 )CALL ELEM08 (R# S, T, SHAPE) 

C 

IF (Ml. EQ. 1 2 ) CALL SHAP12(R, S# T# NCONEC. SHAP, Nl# N3, N4, IJK) 
IF (Ml. EQ. 1 2 )CALL ELEM1 2(R, S. T. SHAPE) 

C 

IF (Ml. EQ. 20 )CALL SHAP20(R, S, T. NCONEC. SHAP, N1.N3.N4, IJK) 
IF (Ml. EQ. 20 )CALL ELEM20(R, S, T, SHAPE) 

C 

DO 20 11-1,3 
DO 20 JJ*1, 3 

SJ( 1 1, JJ)=0. DO 
DO 20 KK=1 , Ml 

I A=NCONEC ( I JK, KK+1 ) 

SJ( 1 1, JJ)=SJ(II, JJ)+SHAPE(II, KK)*COOR(IA, JJ) 

20 CONTINUE 

C####################### ######## 

CALL MULTI P (TRAN, SJ) 

C############################### 

C 

C 

DET=(SJ(2, 2)*SJ(3, 3)-SJ(3, 2)*SJ(2. 3) ) 

C 

c 

DO 30 1=1. INODE 
NIJ=NBOUN( IB, I) 

C 

ELOAD ( I )=ELOAD( I )+SHAP ( NIJ) *D£T*W( J)*W( K) /4. DO 
DO 30 L=1 , INODE 
NIK=NBOUN( IB, L) 

C 

EST IF ( I , L) =ESTIF( I , L >+SHAP ( Nl J) *SHAP (NIK )*D£T* 

> W(J)*W(K)/4. DO 

30 CONTINUE 
C 


40 CONTINUE 


50 



50 CONTINUE 


C 

c 


RETURN 

END 

C 

SUBROUTINE BOUNDS(S/ COOR, NBOUN, NCONEC. ELOAD, IJK, INODE< IB# 
> W, Z,N1,N2,N3,N4, SHAP.ESTIF) 

C 

IMPLICIT REAL*8(A-H» 0-Z > 

DIMENSION COOR (Nl, 3), W(6), Z(6) » SHAPE(3, 20) 

DIMENSION SJC3# 3)» SHAP (Nl ) » NB0UNCN2. 8) » NCONEC (N4» N3+1 ) 
DIMENSION ELOAD (8) # ESTIF(8# 8) 

DIMENSION TRAN <3. 3) 

C 

M1=NC0NEC( IJK, 1 ) 

DO 10 1=1. 8 
ELOAD ( I )=0. DO 
DO 10 1 1=1 #8 
ESTIF( I, II )=0. DO 
10 CONTINUE 
C 

CALL TRANS ( IJK# Nl. N3» N4# COOR# NCONEC. TRAN) 

C 

c 

c 

c 

DO. 50 J= 1.6 

R=0. 5D0*( 1 ; DO+Z ( J) ) 

DO 40 K= 1.6 

T=0. 5D0*C1. DO+Z(K) ) 

C 

C 

IF ( Ml . EQ. 8 ) CALL SHAP08(R. S. T. NCONEC. SHAP. N1.N3. N4. IJK) 

IF (Ml. EQ. 8) CALL ELEM08 (R. S. T. SHAPE) 

C 

IF (Ml. EQ. 12)CALL SHAP 1 2 <R, S, T. NCONEC. SHAP, Nl, N3, N4, IJK ) 
IF (Ml . EQ. 1 2 )CALL ELEM12(R, S, T. SHAPE) 

C 

IF (Ml. EQ. 20)CALL SHAP20(R, S. T. NCONEC. SHAP, Nl, N3. N4, IJK) 
IF (Ml. EQ. 20)CALL ELEM20(R, S, T. SHAPE) 

C 

DO 20 11=1.3 
DO 20 JJ= 1 . 3 

SJ( 1 1, JJ)=0. DO 
DO 20 KK=1, Ml 

IA=NCONEC ( I JK. KK+1 ) 

S J (II, J J ) =S J (II. JJ)+SHAPE( 1 1 . KK) *COOR ( I A, JJ) 

20 CONTINUE 
C 

C######* ##*#*#*# a#################### 

CALL MULT IP ( TRAN, SJ) 

c 
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c 


DET=-SJ< 1 , 3) »SJ(3, 1)+SJ(1, 1)*SJ(3, 3) 

C 

c 

DO 30 1=1, INODE 
NIJ=NBOUN(IB, I) 

C 

ELO AD < I ) -ELOAD ( I ) +SHAP ( NI J ) *D£T*W ( J )*W(K)/4. 00 
DO 30 L«l, INODE 
NIK=NBOUN( IB, L) 

C 

ESTIF(I, L)«ESTIF(I, L)+SHAP(NIJ> *SHAP(NIK)*D£T* 
> W(J)*W(K)/4. DO 

30 CONTINUE 
C 

40 CONTINUE 
50 CONTINUE 
C 

c 

c 


RETURN 

END 

c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
c 


SUBROUTINE BOUNDT ( T, COOR, NBOUN, NCONEC, ELOAD, IJK, INODE, IB. 
> W, Z, N1.N2.N3.N4, SHAP. ESTIF) 


IMPLICIT REAL*8(A-H. 0-Z) 

DIMENSION COOR (Nl , 3), W(6>. Z (6) , SHAPE(3, 20) 

DIMENSION SJ(3. 3), SHAP (Nl >. NBOUN (N2, 8) , NCONEC (N4, N3+1 ) 
DIMENSION ELOAD (8) , ESTIF(8, 8) 

C ################################ 

DIMENSION TRAN (3. 3) 

C #############•## ################# 

c 

M1=NC0NEC( I JK, 1 ) 

DO 10 1=1, 8 
ELOAD ( I )=0. DO 
DO 10 11=1,8 
ESTIF ( I, II )=0. DO 
10 CONTINUE 
C 

C ###*#######################*##**#######«#####« 

CALL TRANS ( IJK. Nl , N3, N4, COOR, NCONEC, TRAN) 

C 

C 

c 

DO 50 J= 1,6 

R=0. 5D0*( 1 . DO+Z ( J) ) 

DO 40 K= 1,6 

S=0. 5D0*( 1 . DO+Z (K) ) 

C 

c 

IF ( Ml . EQ. 8 ) CALL SHAP08(R. S, T, NCONEC. SHAP. N1.N3, N4. IJK) 
IF (Ml. EQ. 8) CALL ELEM08 (R, S. T, SHAPE) 

C 

IF (Ml. EQ. 1 2 )CALL SHAP1 2(R. S. T, NCONEC. SHAP. Nl, N3, N4, IJK > 
IF (Ml . EQ. 12)CALL ELEM12(R. S, T, SHAPE) 

C 
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IF (Ml. EQ. 20)CALL SHAP20(R, S, T, NCONEC, SNAP. N1.N3.N4, I JK ) 
IF (Ml. EQ. 20)CALL ELEM20(R# S. T, SHAPE) 

C 

DO 20 11=1.3 
DO 20 JJ=1 1 3 

SJ(II. JJ)=0. DO 
DO 20 KK=1 . Ml 

I A=NCONEC ( I JK. KK+1) 

SJdl. JJ)=SJ(II. JJ)+SHAPE( 1 1 . KK) *COOR ( IA» JJ) 

20 CONTINUE 
C 

C############################################## 

CALL MULT IP ( TRAN, SJ> 

C 

DET=SJ( 1. 1 )*SJ(2. 2)-SJ(l, 2)*SJ(2. 1) 


DO 30 1=1. INODE 
NIJ=NBOUN(IB, I) 

ELOAD( I )=EL0AD(I)+SHAP(NIJ)*DET*W(J)*W(K>/4. DO 
DO 30 L=1 , INODE 
NIK=NBOUN( IB, L) 

EST IF ( I. L) =ESTIF ( I, D+SHAP ( NIJ) *SHAP (NIK)*D£T* 
> W(J)*W(K)/4. DO 

30 CONTINUE 

40 CONTINUE 
50 CONTINUE 


RETURN 

END 


7.7.*/. X7.7.7. 7.X7.7. X7.7.X 7.7.X7. 7.7.X7.7.7.7.7. 


7. 7.7.7. X7.7.7.7.7.7.7. 7.7.7.7. 7.7. 


SUBROUTINE DECOMP(N, A, LU, INDEXR, NN) 

This routine implements the Gaussian forward elimination 
algorithm to find the LU decomposition of the N by N matrix 
A. Partial pivoting is used along with scaling for row 
equilibration 


IMPLICIT REAL*8 (A-H. O-Z) 

REAL*8 LU 

DIMENSION A(NN, NN) . LU ( NN. NN). SCALEF(200). INDEXR(NN) 

PTOL=. 00000000000000000000000000001 

DO 5 1 = 1, N 

INDEXR ( I )=I 
SCALEF ( I )=A (I, 1 ) 

DO 4 J= 1 , N 

LU( I » J)=A( I, J) 

IF ( DABS (LU( I. J) ) . GT. SCAuEF(I)) SC ALEF (I >=DABS ( LU ( I , J) ) 
4 CONTINUE 
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c 

c 

IF (SCALEF(I). NE. 0. DO) GO TO 9 
WRITE <2, *> 'ALL ZERO ROW' 

RETURN 
5 CONTINUE 
C 

c 

NM 1=N- 1 
DO 50 K-l. NM1 

BIG=DABS(LU (K. K)/SCALEF(K) ) 

IBIG=K 

KP1=K+1 

DO 10 IR=KP1,N 

IF(DABS( LU( IR. K) /SCALEF( IR) ) . LE. BIG) GOTO 10 
BIG=DABS(LU<IR. K)/8CALEF<IR> > 

IBIG-IR 

10 CONTINUE 

C 
C 

IF(BIG. GT. PTOL) GOTO 12 
WRITE(2,*) 'SMALL PIVOT' 

RETURN 

12 IF( IBIG. EQ. K) GOTO 16 

C 
C 

ISAVE»INDEXR<K> 

INDEXR (K)=INDEXR(IBIG) 

INDEXR ( IBIG )®ISAVE 
SAVE=SCALEF (K) 

SCALEF ( K ) ss SCALEF (IBIG) 

SCALEF (IBIG )=SAVE 
C 

c 

DO 15 J-l. N 

SAVE=LU< IBIG. J) 

LU( IBIG. J)=LU<K. J) 

LU( K. J)=SAVE 

15 CONTINUE 
C 

c 

16 DO 30 I =KP 1 . N 

C 

IF(LU(I. K). EQ. 0. DO) GO TO 30 
LU< I • K)**LU( I . K) /LU(K. K) 

DO 20 J*KP1.N 

LU(I. J)*LU(I. J)-LU(I. K)*LU(K. J) 

20 CONTINUE 

30 CONTINUE 

50 CONTINUE 
C 

c 

IF <DABS(LU (N. N ) /SCALEF (N) ) . GT. PTOL) GO TO 60 
WRITE (2,*) 'SMALL PIVOT' 

60 RETURN 
END 
C 
C 


c y.v.v.v, y.'/:/:/. v.y.y.y. y.y.y.y. y.y.y.y. y.y.y.y.y.y.%y. y.xy.y. y.y.y.y. y.y.y.y. y.y.y.y. */.•/.■/.•/. 


c 
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C 

SUBROUTINE SOLVE(N, LU, B. X, INDEXR, NN) 

This routine solves the linear equations with right-hand 
side vector B using the LU decomposition matrices already 
computed by routine decomp. The final solution is stored 
in X array. 


IMPLICIT REAL*B (A-H, O-Z) 

REAL*8 LU 

DIMENSION LU!NN, NN), B(NN), X(NN), INDEXR(NN) 


DO 1 1 = 1, N 

X!I )=B! INDEXR(I ) ) 
1 CONTINUE 


NM1=N-1 
DO 3 1=2, N 
IM1 =1-1 
DO 2 J=l, IM1 

X<I)=X(I ) -LU ( I, J )*X ( J) 

2 CONTINUE 

3 CONTINUE 


X(N)=X(N)/LU(N, N) 

DO 5 11 = 1, NM1 
I=N-I I 
IP1=I+1 
DO 4 J=IP1, N 

X < I > =X < I )-LU( I, J)*X( J) 

4 CONTINUE 

X(I )=X ( I ) /LU < I » I) 

5 CONTINUE 
RETURN 
END 


7 . 7 . v :/:/:/. xtosx v.v.y.7.v.7.y.v. 7.v.7.v.v.7.7.v. xxxx v.v.v.v . xxxx */.*/. 


SUBROUTINE ELEM08(R, S, T. SHAPE) 
IMPLICIT REAL* 8 (A-H, O-Z) 
DIMENSION SHAPE (3, 20) 

SHAPE! 1, 1 ) = <1. DO-S )*( 1 . DO-T) 
SHAPE (1,2) =S* < 1 . DO-T) 

SHAPE (1,3) =S* ( T-l . DO) 

SHAPE ( 1, 4) = <S-1. D0)*( 1. CO-T) 
SHAPE! 1, 5) = ( 1. DO-S ) *T 
SHAPE! 1,6) =S*T 
SHAPE (1,7) =— 1 . DO*S*T 
SHAPE! 1, 8) = (S-1. DO ) *T 


SHAPE (2,1) =R* ( T-l . DO) 
SHAPE (2,2) =R* ( 1. DO-T) 
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c 

c 


c 


c 

c 

c 

c 


c 


c 

c 


c 

c 


SHAPE ( 2, 3) = ( 1. DO-R ) * ( 1 . DO-T) 
SHAPE (2,4)* (R- 1 . DO ) * ( 1 . DO-T) 
SHAPE ( 2, 5) =-l. DO*R*T 
SHAPE (2, 6) =R*T 
SHAPE ( 2i 7) = ( 1. DO-R ) *T 
SHAPE ( 2/ B ) = (R— 1 . DO )*T 


SHAPE (3, 1 ) =R*( S-l. DO) 

SHAPE ( 3, 2) «-l. DO*R*S 
SHAPE < 3, 3)»(R-1. DO)*S 
SHAPE ( 3, 4) = (R-1. DO )*( 1 . DO-S) 
SHAPE ( 3f 9 ) «R* ( 1 . DO-S) 

SHAPE (3, 6)*R#S 

SHAPE ( 3, 7) = ( 1 . DO-R ) *S 

SHAPE ( 3, 8) *( 1. DO-R )*( 1 . DO-S) 

RETURN 

END 


%y. 7.V.Y.7. v.v.y.v. v.y.%y. y.w.v. ■///.%% v.v.v.y. v.y.%7. xxxx y.yy.v. xxxx xxxx xxxx xxxx 


SUBROUTINE ELEM12( R, S. T, SHAPE) 
IMPLICIT REAL*8 (A-H, O-Z) 
DIMENSION SHAPE(3, 20) 


SHAPE ( 1, 1> = (1. DO-S )*< 1 . DO-T) 

SHAPE ( 1, 2>=S*<2. D0*T*T-4. D0*R*T+4. D0*R-T-1. DO) 

SHAPE ( 1, 3)=S*(-2. D0*T*T-4. D0*R*T+4. D0*R+5. D0*T-3. DO) 
SHAPE ( 1, 4) = (S-1. D0)*(1. DO-T) 

SHAPE ( 1/ 5) = ( 1. DO-S )*T 

SHAPE < 1 , 6 ) =S* ( 2. D0*T*T+4. D0*R*T-3. DO*T ) 

SHAPE (1,7) =S* ( 4. D0*R*T-2. DO*T*T-T ) 

SHAPE ( 1, 8)=(S-1. DO)*T 

SHAPE ( 1, 9> = <4. DO-8. D0*R)*S*(1. DO-T) 

SHAPE ( 1, 10 )=S*4. DO*T*( 1. DO-T) 

SHAPE ( 1, 11 )=S*4. DO*T*( T-l. DO) 

SHAPE ( 1, 12 ) = (4. DO-B. DO*R)*S*T 


SHAPE (2, 1 ) =R* ( T-l. DO) 

SHAPE (2, 2) =2. DO*R*T*T— 2. D0*R*R*T+2. DO*R*R-R*T-R 
SHAPE (2, 3) =2. D0*(R*R+T*T-R*T*T-R*R*T)+3. DO*R*T 
> -3. D0*T-3. D0*R+1. DO 

SHAPE (2, 4 ) = (R— 1. D0)*(1. DO-T) 

SHAPE (2, 5) =-l. DO*R*T 

SHAPE (2, 6 ) - 2 . D0*R*T*T+2. D0*R*R*T-3. DO*R*T 
SHAPE (2, 7) =2. D0*R*R*T-2. D0*R*T*T+2. DO*T*T-R*T-T 
SHAPE < 2, 8) = (R-1. DO )*T 
SHAPE ( 2, 9) =4. D0*R*<1. D0-R)*(1. DO-T) 

SHAPE (2, 10 )-R*4. DO*T*< 1. DO-T) 

SHAPE ( 2, 11 )-(l. DO-R) *4. D0*T*<1. DO-T) 

SHAPE <2, 12) =4. D0*R*<1. DO-R)*T 


SHAPE (3, 1 ) =R* < S-l . DO) 

SHAPE ( 3, 2 ) =S* ( 4. D0*R*T-2. DO*R*R-R ) 

SHAPE (3,3) =S* < -4. D0*R*T-2. D0*R*R+4. DO*T+S. D0*R-3. DO) 
SHAPE (3, 4 ) “ (R — 1. D0)*(1. DO-S) 
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c 

c 

c 


c 


c 

c 


c 


SHAPE (3,5) =R* ( 1 . DO-S) 

SHAPE (3,6) =S* ( 4. DO*R*T +2. D0*R*R-3. DO*R ) 

SHAPE (3,7) =S* ( 2. D0*R*R-4. D0*R*T+4. D0*T-R-1. DO) 
SHAPE (3, 8) = ( 1. DO-R > * ( 1 . DO-S) 

SHAPE ( 3, 9) =4. DO*R* <R— 1 . DO)*S 
SHAPE (3, 10 )=R*S*(4. DO-8. DO*T) 

SHAPE ( 3, 11 ) = < 1 . DO-R )*S*(4. DO-8. DO*T) 

SHAPE (3, 12) =4. D0*R*( 1. DO-R)*S 


RETURN 

END 


y.7.x7.7.7.y.7.’///. xxxxxxxx y.7.7.%7.7.y.xz%%%xx?:xx%xxy.%xxx7.xy.7.y.y.y.*/.’/.5iy.%y. 


SUBROUTINE ELEM20( R, S, T, SHAPE) 

IMPLICIT REAL* 8 (A-H, 0-Z) 

DIMENSION SHAPE(3, 20) 

SHAPE ( 1, 1 >e(l-S)*( 1-T)*(4*R-2*S-2*T-1> 
SHAPE ( 1, 2)=S*( 1-T) *(4*R+2*S-2*T-3) 
SHAPE ( 1, 3)=S*( 1-T)*(4*R-2*S+2*T-1) 
SHAPE ( 1, 4) = <1-S)*( 1-T) *<4*R+2*5+2*T-3) 
SHAPE ( 1 , 5 ) = ( 1 -S ) *T* ( 4*R-2*S+2*T-3 ) 
SHAPE ( 1, 6) =>S*T*(4*R+2*S+2*T-5) 

SHAPE ( 1, 7)=S*T*(4*R-2*S-2*T+1 ) 

SHAPE ( 1, 8) = (1-S)*T*(4*R+2*S-2*T-1 ) 
SHAPE (1,9) =4*S* ( 1— S)* ( 1-T) 

SHAPE ( 1, 10 )=<4-8*R )*S* ( 1-T ) 

SHAPE ( 1 , 11 )=4*S*(S-1)*(1-T> 

SHAPE ( 1 , 12>=(4-8*R)*(1-S)*(1-T) 

SHAPE < 1, 13 )=4*S*( 1-S)*T 
SHAPE ( 1 , 1 4 ) = ( 4-8*R ) *S*T 
SHAPE ( 1, 15) =4*S*( S— 1 ) *T 
SHAPE ( 1, 16 ) = (4-B*R )*(1-S)*T 
SHAPE ( 1, 17 ) =( 1 -S) *4*T* ( 1— T ) 

SHAPE ( 1, 18 ) =S*4*T* (1-T ) 

SHAPE ( 1, 19) =S*4*T*(T-1 ) 

SHAPE ( 1 , 20 ) = (S-1 )*4*T* ( 1-T ) 


SHAPE (2,1) =R* ( 1 -T ) * ( -2*R+4*S+2*T— 1 > 
SHAPE (2, 2) =R*( 1-T) * (2*R+4*S-2*T-3) 

SHAPE ( 2, 3) =< 1-R )*( 1-T)*(-2*R+4*S-2*T-1 ) 
SHAPE (2, 4 ) = ( 1 -R ) * ( 1-T) *(2*R+4*S+2*T-3) 
SHAPE (2, 5)=R*T*<-2*R+4*S-2*T+1) 

SHAPE (2, 6) =R*T*(2*R+4*S+2*T— 5) 

SHAPE (2, 7) = (1-R)*T*( -2*R+4 *S+2*T-3 ) 
SHAPE ( 2, 8 ) = ( 1 -R ) *T* (2*R+4*S-2*T— 1 ) 

SHAPE (2, 9>=R*(4-8*S)*< 1-T) 

SHAPE ( 2, 10)=4*R*(1-R)*(1-T> 

SHAPE (2, 11 )= S (1-R)*( 4-8 *S) * ( 1-T ) 

SHAPE ( 2, 12 )=4*R*(R— 1 )* ( 1-T ) 

SHAPE (2, 13 ) -R* ( 4-B*S) *T 
SHAPE ( 2, 14)=4*R*(1-R)*T 
SHAPE ( 2, 15 ) = ( 1-R)* (4-B*S)*T 
SHAPE ( 2, 16)=4*R*(R-1 )*T 
SHAPE (2, 17)=R*4*T*(T-1 ) 

SHAPE <2, 18 ) =R*4*T* ( 1-T ) 

SHAPE ( 2, 19 ) = ( 1 -R ) *4*T * ( 1-T ) 

SHAPE (2, 20 ) = < R-l ) *4*T* ( 1— T ) 
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SHAPE ( 3* 1 )»R*< 1-S) *(— 2*R+2*S+4*T-1 ) 
SHAPE ( 3» 2 ) =R*S*<-2*R-2*S+4*T+1 ) 

SHAPE < 3. 3 ) * ( 1 -R ) *S* (2*R-2*S+4*T-1 ) 
SHAPE ( 3, 4 ) = ( 1— R )* ( 1— S) * (2*R+2*S+4*T— 3) 
SHAPE (3. 5>»R*< l-S)*<2*R-2*S+4*T-3) 

SHAPE (3. 6) =R*S*(2*R+2*S+4*T-5) 

SHAPE ( 3, 7 ) = ( 1-R ) *S*(-2*R+2 *S+4 *T-3 ) 
SHAPE (3, 8 ) = ( 1-R )* ( 1-S) *(-2*R-2*S+4*T— 1 ) 
SHAPE (3, 9)=R*4*S*(S-1 ) 

SHAPE ( 3, 10)=4*R*(R-1)*S 
SHAPE ( 3, 11 )=<R-1)*4*S*(1-S) 

SHAPE (3, 12)=4*R*(R-1)*(1-S) 

SHAPE (3,13) =R*4*S* ( 1-S ) 

SHAPE (3, 14) =4*R* ( 1 -R ) *S 
SHAPE ( 3, 1 5 ) = ( 1-R) *4*S* < 1-S ) 

SHAPE ( 3, 16 )=4*R*( 1-R)* ( 1-S ) 

SHAPE < 3* 17 )=R* < 1-S )*(4-8*T ) 

SHAPE (3, 18 )=R*S*(4-B*T ) 

SHAPE ( 3, 19 )=< 1 -R)*S*(4-8*T ) 

SHAPE (3, 20 >=( 1 -R ) * ( 1-S ) *(4-8*T ) 

RETURN 

END 


7.7. 7. y.7. 7. 7.7.7 7. 77.7.7. 7.7.77 7777 7.7.77 77.77. 7777 7 7.7.7 777.7 7.7.7// 7.77.7. 7.7.77. 7.7.77 X777. 


SUBROUTINE SHAP08 ( R, S» T# NCONEC, SHAP, N1 , N3, N4, IJK) 

IMPLICIT REAL* 8 (A-H, 0-Z) 

DIMENSION NC0NEC(N4, N3+1), SHAP(Nl) 


SHAP (NCONEC ( IJK, 2 ) ) =R* ( 1-S ) *( 1 -T> 
SHAP < NCONEC ( IJK, 3) )=R*S*(1-T) 

SHAP (NCONEC (IJK, 4) )=( 1-R)*S*( 1-T) 

SHAP (NCONEC ( IJK, 5 ) )=( 1 -R) * ( 1-S )*( 1 -T) 
SHAP (NCONEC (IJK, 6) )=R*(1-S)*T 
SHAP ( NCONEC ( I JK, 7 ) ) =R*S*T 
SHAP (NCONEC (IJK, 8) )=(1-R)*S*T 
SHAP ( NCONEC ( I JK, 9 ) )=( 1 -R) * ( 1-S )*T 

RETURN 

END 


C 


C 


XX 7 7.7.7. 7.7. 77 7.7.7.X 7X7.7. 7777 7777 7777 7777 7777 7777 7777 7777 y m y. m / m y. m /.7. m / m 


SUBROUTINE SHAP12(R,S, T, NCOhEC, SHAP, Nl, N3, N4, IJK) 

IMPLICIT REAL*B (A-H, 0-Z) 

DIMENSION NCONEC <N4, N3 + 1), SHAP(Nl) 


SHAP (NCONEC (IJK, 2) )=R*(1-S)*< 1-T) 

SHAP ( NCONEC ( I JK, 3 ) ) =S* ( 2*R *T*T -2*R *R*T +2*R *R-R * T-R ) 

SHAP ( NCONEC (IJK, 4) ) =S* ( -2*R*T*T-2*R*R*T+2*R*R+2*T*T+5*R*T 
> -3*T-3*R+1 ) 

SHAP < NCONEC ( I JK, 5 ) ) = ( 1 -R ) * ( 1-S ) *( 1 -T) 

SHAP (NCONEC (IJK, 6) )=R*( 1-S)*T 

SHAP (NCONEC (IJK, 7) )=S* <2*R*T*T+2*R*R*T-3*R*T) 
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c 


SHAP ( NCONEC ( I JK, 8 ) > =S* ( 2*R *R*T-2*R*T*T +2*T *T-R * T-T ) 
SH AP ( NCONEC ( I JK, 9 ) ) = ( 1-R)*(1-S)*T 
SHAP ( NCONEC < I JK, 10) )=4*R*( 1-R)*S*< 1-T) 

SHAP (NCONEC (IJK, 11 ) )=R*S*4*T*< 1-T) 

SHAP ( NCONEC ( IJK, 12) )=< 1-R)*S*4#T*< 1-T) 

SHAP (NCONEC (IJK, 13) )=4*R*( 1-R)*S*T 


RETURN 

END 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX9 IXXXXXXXXXXXXXXXXXXXXX 


SUBROUTINE SHAP20(R, S, T, NCONEC# SHAP.N1# N3# N4» I JK) 

IMPLICIT REAL*8 <A-H#0-Z) 

DIMENSION NCONEC (N4» N3+1 ) # SHAP (N1 ) 


SHAP ( NCONEC ( I JK# 2) )=R* ( 1-S )*( 1-T)* <2*R-2*S-2*T-1 ) 

SHAP < NCONEC ( I JK, 3) )=R*S*( 1 -T) * <2*R+2*S-2*T-3) 

SHAP ( NCONEC ( I JK, 4 ) ) = ( 1 -R ) *S*< 1 -T) * <-2*R+2*S-2* T-l ) 
SHAP ( NCONEC ( I JK, 5 ) ) = ( 1 -R ) * ( 1-S ) *< 1 -T) * ( 1-2*R-2*S-2*T ) 
SHAP ( NCONEC ( I JK, 6 ) )=R* ( 1-S )*T* (2*R-2*S+2*T-3) 

SHAP (NCONEC (IJK, 7) )=R*S*T* <2*R+2*S+2*T-5) 

SHAP ( NCONEC (IJK,8)) = (1-R) *S*T* < -2*R+2*S+2*T-3 ) 

SHAP (NCONEC ( I JK, 9 ) ) = ( 1 -R ) * < 1-S ) *T* ( -2*R-2*S+2* T-l ) 
SHAP (NCONEC (IJK, 10 ) )=R*4*S*( 1-S)*( 1-T) 

SHAP (NCONEC (IJK, 11 ) )=4*R*( 1-R)*S*< 1-T) 

SHAP (NCONEC (IJK, 12) )=( 1-R ) *4*S*( 1-S)*< 1-T) 

SHAP (NCONEC (IJK, 13) )=4*R*( 1-R)*(1-S)*< 1-T) 

SHAP (NCONEC (IJK, 14 ) )=R*4*S*( 1-S)*T 
SHAP ( NCONEC ( IJK, 15) )=4*R*( 1-R)*S*T 
SHAP (NCONEC (IJK, 16) )=( 1-R ) *4*S*( 1-S)*T 
SHAP (NCONEC (IJK, 17) )=4*R*( 1-R)*(1-S)*T 
SHAP ( NCONEC ( I JK, 18) )=R*( 1-S)*4*T*( 1-T) 

SHAP (NCONEC (IJK, 19) )=R*S*4*T*( 1-T) 

SHAP(NCONEC< IJK, 20) )=( 1-R) *S*4*T*( 1-T) 

SHAP (NCONEC (IJK, 21 ) >=( 1-R ) *( 1-S)*4*T*< 1-T) 


RETURN 

END 




SUBROUTINE MULT IP ( TRAN. SJ) 

IMPLICIT REAL*8(A-H, 0-Z) 

DIMENSION TRAN (3, 3 ) , SJ ( 3, 3 ) , CH ( 3, 3 ) 

DO 10 1 = 1,3 
DO 10 J=l, 3 

CH( I, J)=0. DO 
DO 10 K=1 , 3 

CH( I, J) =CH( I, J)+TRAN(I. K)*SJ(K# J) 
10 CONTINUE 

DO 20 1 = 1, 3 
DO 20 J=l, 3 

SJ ( I , J)=CH( I, J) 

20 CONTINUE 
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RETURN 

END 


SUBROUTINE TRANSd.Nl, N3, N4, COOR, NCONEC. TRAN) 
IMPLICIT REAL*B(A-H» O-Z) 

DIMENSION COOR (Nl» 3), NCONEC <N4, N3+1 ), TRAN (3, 3) 


XO=COOR (NCONEC ( I, 5 ) . 1 ) 
YO=COOR (NCONEC < I, 5 ) , 2 ) 
ZO=COOR (NCONEC (I, 5), 3) 

C 

X 1 =COOR ( NC ONEC (1.2). l)-XO 
Y1 =COOR (NCONEC (1.2) , 2) -YO 
Z 1 “COOR ( NC ONEC ( 1 . 2),3)-Z0 
C 


X2=C00R (NCONEC ( I, 4 ) . 1 ) -XO 
Y2=C00R ( NC ONEC ( I. 4 ). 2) -YO 
Z2=C00R (NCONEC (1.4). 3) -ZO 
C 

X3=C00R (NCONEC (1,9), l)-XO 
Y3=C00R (NCONEC ( I. 9 ) . 2) -YO 
Z3=C00R ( NC ONEC ( I, 9), 3) -ZO 
C 

A1=DSQRT ( X 1*X 1 +Y1*Y1+Z 1*Z1 ) 
A2=DSQRT( X2*X2+Y2*Y2+Z2*Z2 ) 
A3=DSQRT(X3*X3+Y3*Y3+Z3*Z3) 
C 

TRAN( 1, 1 >=X1/A1 
TRAN< 1, 2) = X2/A2 
TRAN(1, 3) = X3/A3 
C 

TRAN<2, 1 )=Y1/A1 
TRAN(2. 2) = Y2/A2 
TRAN(2. 3 ) =Y3/A3 
C 

TRAN( 3. 1 >=Z1/A1 
TRAN( 3. 2>=Z2/A2 
TRAN( 3. 3) =Z3/A3 
C 


RETURN 

END 


C 

C 


c •/.%*/. v.y.7.7. 'a;/.'/. 7.7.y.7.y.v.v.y. 7 . 7 . 7 .% xxxx v.v.y.7. v.v.v.7. y.7.xv. y.v.y.7. xxxx xxxx %'/:/'/. xx%xxxxx v.y.7.% v. 
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